# Searching for Lindbladians obeying local conservation laws and showing
thermalization

Devashish Tupkary<EMAIL_ADDRESS>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>Department of Physics, Indian Institute of Technology, Hyderabad 502285, India
Centre for complex quantum systems, Aarhus University, Nordre Ringgade 1, 8000
Aarhus C, Denmark  School of Physics, Trinity College Dublin, College Green,
Dublin 2, Ireland

###### Abstract

We investigate the possibility of a Markovian quantum master equation (QME)
that consistently describes a finite-dimensional system, a part of which is
weakly coupled to a thermal bath. In order to preserve complete positivity and
trace, such a QME must be of Lindblad form. For physical consistency, it
should additionally preserve local conservation laws and be able to show
thermalization. First, we show that the microscopically derived Redfield
equation (RE) violates complete positivity unless in extremely special cases.
We then prove that imposing complete positivity and demanding preservation of
local conservation laws enforces the Lindblad operators and the lamb-shift
Hamiltonian to be ‘local’, i.e, to be supported only on the part of the system
directly coupled to the bath. We then cast the problem of finding ‘local’
Lindblad QME which can show thermalization into a semidefinite program (SDP).
We call this the thermalization optimization problem (TOP). For given system
parameters and temperature, the solution of the TOP conclusively shows whether
the desired type of QME is possible up to a given precision. Whenever
possible, it also outputs a form for such a QME. For a XXZ chain of few
qubits, fixing a reasonably high precision, we find that such a QME is
impossible over a considerably wide parameter regime when only the first qubit
is coupled to the bath. Remarkably, we find that when the first two qubits are
attached to the bath, such a QME becomes possible over much of the same
paramater regime, including a wide range of temperatures.

##  I Introduction

A small finite-dimensional quantum system, a part of which is weakly coupled
to a macroscopic thermal bath, is expected to thermalize to the temperature of
the bath. Describing this dynamics is relevant across various fields in
quantum science and technology, including quantum information and
thermodynamics Binder  _et al._ (2018), quantum optics Hur  _et al._ (2016),
quantum chemistry Cao  _et al._ (2019), engineering Awschalom  _et al._ (2021)
and biology Cao  _et al._ (2020). In absence of coupling to the macroscopic
thermal bath, the dynamics of the density matrix of the system is governed by
the Heisenberg equation of motion. This unitary evolution is Markovian. When
coupled to the macroscopic thermal bath, the dynamics becomes non-unitary,
described by a quantum master equation (QME) Breuer and Petruccione (2006);
Rivas and Huelga (2012); Carmichael (2002). We investigate whether it is
possible to have a physically consistent Markovian QME describing such
dynamics. In order to do so, we are led to introduce the “thermalization
optimization problem” (TOP). This is a semidefinite program (SDP), the output
of which conclusively shows whether, for given system parameters and
temperature, such a QME is possible, up to given precision. Whenever possible,
the output also yields one possible form for such a QME. Whenever impossible,
it means that, for such parameters, the dynamics cannot be described by any
Markovian QME even at weak system-bath coupling, and therefore must have some
non-Markovian character. The SDP can be solved using standard packages in
high-level computing. We note that, while SDP is widely used in many branches
of quantum information and communication Watrous (2018); Siddhu and Tayur
(2022), and also in quantum chemistry Mazziotti (2004, 2011), it has been
combined with open quantum system techniques in only few previous works Wolf
_et al._ (2008); Cubitt  _et al._ (2012); Kiilerich and Mølmer (2018), in very
different contexts.

Figure 1: Schematic of the setup we consider: an arbitrary finite dimensional
system described by Hamiltonian $H_{S}$, a part of which is weakly coupled to
a thermal bath at inverse temperature $\beta$. The Hilbert space of the
system, $\mathcal{H}_{S}$ is divided into a part $\mathcal{H}_{L}$ which
directly couples to the bath, and the remaining part $\mathcal{H}_{M}$.

It was shown by Gorini, Kossakowski, Sudarshan, and Lindblad (GKSL) Gorini
_et al._ (1976, 1978); Lindblad (1976) that any QME that preserves complete
positivity and trace of the density matrix, and describes Markovian dynamics
has to be of the form

$\displaystyle\frac{\partial\rho}{\partial
t}=i[\rho,H_{S}+H_{LS}]+\mathcal{D}(\rho),$ (1)
$\displaystyle\mathcal{D}(\rho)=\sum_{\lambda=1}^{d^{2}-1}\gamma_{\lambda}\Big{(}L_{\lambda}\rho
L_{\lambda}^{\dagger}-\frac{1}{2}\\{L_{\lambda}^{\dagger}L_{\lambda},\rho\\}\Big{)},~{}~{}\gamma_{\lambda}\geq
0,$

which is commonly called the Lindblad equation. In Eq. (LABEL:Lindblad_form),
$\rho$ is the density matrix of the system, $d$ is the Hilbert space dimension
of the system, $H_{S}$ is the system Hamiltonian, $H_{LS}$ is the Lamb shift
Hamiltonian, $L_{\lambda}$ are the Lindblad operators, $\gamma_{\lambda}$ are
the rates, and $\mathcal{D}(\rho)$ is called the “dissipator” term. The
preservation of complete positivity condition is enforced by demanding
$\gamma_{\lambda}\geq 0$. Lindblad equations have been extensively used in
studying both theoretical and experimental setups Breuer and Petruccione
(2006); Carmichael (2002); Rivas and Huelga (2012); Weimer  _et al._ (2021);
Plenio and Knight (1998); Dalibard  _et al._ (1992); Mølmer  _et al._ (1993).

Given this enormous scope of application, it is of paramount importance to
assess the conditions under which such a Markovian description emerges from a
more microscopic theory. The standard way to microscopically obtain a
Markovian QME is to consider the global Hamiltonian of the system weakly
coupled to baths, and to trace out the baths perturbatively to the leading
order. Starting with this microscopic viewpoint, it becomes clear that only
having an equation in the Lindblad form is not sufficient, there are some
additional fundamental requirements for physical consistency Tupkary  _et al._
(2022). In particular, one must preserve local conservation laws, and if the
steady state is unique, the system is not driven and all baths have same
temperature, the system must thermalize to the temperature of the baths. It
would be useful to have a QME which, by construction, is of Lindblad form and
satisfy these additional requirements. In this paper, we systematically go
about searching for such a QME for a setup where a part of the system is
coupled to a single bath (see Fig. 1). This is done in three steps, each step
having important consequences:

  1. 1.

The microscopically derived quantum master equation to the leading order in
system-bath coupling is the so-called Redfield equation (RE) Breuer and
Petruccione (2006). The RE has been shown to preserve local conservation laws
and be able to show thermalization Tupkary  _et al._ (2022). Here, we provide
an explicit, model independent proof that the RE necessarily violates complete
positivity unless the Redfield dissipator happens to act “locally”, meaning it
is identity on the part of the system that is not directly coupled to the
bath. Although the violation of complete positivity by RE has been previously
demonstrated in specific examples 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), we are
unaware of any previous work with such a model independent explicit proof.

  2. 2.

We then prove that enforcing complete positivity condition
$\gamma_{\lambda}\geq 0$ and preservation of local conservation laws
necessarily requires the Lindblad operators and the Lamb shift Hamiltonian to
be “local”. That is, they must act only on the part of the system coupled to
the bath, and be identity on the part of the system that is not connected to
the bath. This directly rules out the possibility of any so-called ‘global’
Lindblad equation, such as the eigenbasis Lindblad equation Breuer and
Petruccione (2006); Carmichael (2002), the Universal Lindblad Equation Nathan
and Rudner (2020) to be consistent with local conservation laws.

  3. 3.

Given the restriction of the dissipator and the Lamb shift Hamiltonian to be
“local’, we devise a numerical technique using SDP to check conclusively in a
case-by-case basis whether such a QME can show thermalization in a particular
situation. We call this the TOP. We use this method to study the case of a XXZ
chain of few qubits with a part of it coupled to a bath. If the bath is
coupled only to the first qubit, our method conclusively shows that over a
large regime of system parameters and temperature, no such QME exists.
However, remarkably if the bath is coupled to two qubits of the chain, up to a
chosen precision, our method shows that a Marovian QME respecting all
conditions becomes possible over a considerable regime of parameters,
including a wide range of temperatures. Note that the RE for the XXZ chain
does show thermalization and preserve local conservation laws Tupkary  _et
al._ (2022), even when only one qubit is attached to a bath. But it is not
completely positive.

This work is organized as follows. In Sec. II we explain the setup studied in
this work, derive the Redfield equation for our setup, and show that it will
necessarily violate complete positivity, unless the Redfield dissipator
happens to act “locally”. In Sec. III we consider quantum master equations
preserving complete positivity and obeying local conservation laws, and show
that such equations must have a dissipator and Lamb shift operator that acts
only on the part of the system coupled to the bath. In Sec. IV we discuss the
possibility of QMEs respecting complete positivity, local conservation laws
and being able to show thermalization. We introduce the TOP, and use it in the
specific case of the few site XXZ chain with one or two sites attached to
bath. In Sec. V we summarize our results, and discuss future directions.
Certain details are delegated to the appendices.

##  II The Model and Redfield Description

Our setup is described schematically in Fig. 1. The full Hamiltonian be given
by

$H=H_{S}+\epsilon H_{SB}+H_{B},$ (2)

where $H_{S}$ and $H_{B}$ are the Hamiltonians of the system and bath
respectively, $\epsilon\ll 1$ is a small dimensionless parameter controlling
system-bath coupling strength, and $H_{SB}$ is the system-bath coupling
Hamiltonian. At initial time, the system is considered to be in an arbitrary
initial state $\rho(0)$, while the bath is in a thermal state with inverse
temperature $\beta$

$\displaystyle\rho_{\rm
tot}(0)=\rho(0)\otimes\rho_{B},~{}~{}\rho_{B}=\frac{e^{-\beta H_{B}}}{Z_{B}}.$
(3)

Starting with this initial state, the whole set-up of the system and the bath
is evolved with the full Hamiltonian $H$, and the bath degrees of freedom are
traced out to obtain the state of the system,

$\displaystyle\rho(t)={\rm Tr}_{B}\left(e^{-iHt}\rho_{\rm
tot}(0)e^{iHt}\right),$ (4)

where ${\rm Tr}_{B}(\ldots)$ denotes trace over bath degrees of freedom. The
Eq. (4), by construction, is a completely positive trace preserving (CPTP) map
from $\rho(0)$ to $\rho(t)$ Breuer and Petruccione (2006); Rivas and Huelga
(2012). Without any loss of generality, we assume ${\rm
Tr}_{B}(H_{SB}\rho_{B})=0$, where ${\rm Tr}_{B}(\ldots)$ denotes trace over
bath degrees of freedom Breuer and Petruccione (2006); Rivas and Huelga
(2012). The effective equation of motion for the system density matrix written
to the leading order in system-bath coupling strength $\epsilon$ is the RE,
given by, Breuer and Petruccione (2006),

$\displaystyle\frac{\partial\rho}{\partial t}=$ $\displaystyle
i[\rho(t),H_{S}]$ (5) $\displaystyle+$
$\displaystyle\epsilon^{2}\int_{0}^{\infty}dt^{\prime}\operatorname{Tr}_{B}[H_{SB},[H_{SB}(-t^{\prime}),\rho(t)\otimes\rho_{B}]],$

where

$H_{SB}(t)=e^{i(H_{S}+H_{B})t}H_{SB}e^{-i(H_{S}+H_{B})t}$ (6)

and $\rho_{B}$ is the state of the bath. In complete generality, we can write
the system-bath coupling Hamiltonian as

$H_{SB}=\sum_{l}(S_{l}B_{l}^{\dagger}+S_{l}^{\dagger}B_{l}),$ (7)

where $S_{l}$ and $B_{l}$ are operators on the system and bath respectively,
and $l$ can be summed over as many indices as required for $H_{SB}$. Using
Eq.(7) in Eq.(5) and simplifying, we have

$\displaystyle\frac{\partial{\rho}}{\partial t}=$
$\displaystyle\,i[\rho(t),H_{S}]+\epsilon^{2}\Big{\\{}$
$\displaystyle\sum_{l}\big{[}S^{\dagger}_{l},S_{l}^{(1)}\rho(t)\big{]}-\big{[}S^{\dagger}_{l},\rho(t)S_{l}^{(2)}\big{]}$
(8) $\displaystyle+\text{H.c}\Big{\\}},$

where

$\displaystyle S_{l}^{(1)}$
$\displaystyle=\sum_{m}\int_{0}^{\infty}dt^{\prime}\operatorname{Tr}\big{[}B_{l}B^{\dagger}_{m}(-t^{\prime})\rho_{B}\big{]}S_{m}(-t^{\prime})$
(9)
$\displaystyle+\sum_{m}\int_{0}^{\infty}dt^{\prime}\operatorname{Tr}\big{[}B_{l}B_{m}(-t^{\prime})\rho_{B}\big{]}S^{\dagger}_{m}(-t^{\prime})$
$\displaystyle S_{l}^{(2)}$
$\displaystyle=\sum_{m}\int_{0}^{\infty}dt^{\prime}\operatorname{Tr}\big{[}B^{\dagger}_{m}(-t^{\prime})B_{l}\rho_{B}\big{]}S_{m}(-t^{\prime})$
$\displaystyle+\sum_{m}\int_{0}^{\infty}dt^{\prime}\operatorname{Tr}\big{[}B_{m}(-t^{\prime})B_{l}\rho_{B}\big{]}S^{\dagger}_{m}(-t^{\prime}).$

Since the actual microscopic evolution is given by a CPTP map [see Eqs. (3),
(4)], one might naively expect that the evolution obtained from the
microscopically derived RE respects complete positivity. However, as we prove
in the next subsection in generality, unless in extremely special cases, the
RE violates complete positivity.

###  II.1 Violation of complete positivity in Redfield equation

####  II.1.1 Choosing an operator Basis

As shown in Fig. 1, we consider only a part of the system is coupled to the
bath. Let us denote $\mathcal{H}_{L}$ as the Hilbert space of that part of the
system that couples to the bath, and let $\mathcal{H}_{M}$ be the Hilbert
space of the remaining part of the system. In mathematical terms, this means
that any operator $O_{M}$ in the Hilbert space $\mathcal{H}_{M}$ commutes with
the system-bath coupling Hamiltonian $H_{SB}$,

$\displaystyle[O_{M},H_{SB}]=0.$ (10)

The system Hamiltonian can then be written as

$\displaystyle H_{S}=H_{L}+H_{LM}+H_{M},$ (11)

where the Hamiltonian $H_{L}$ is in Hilbert space $\mathcal{H}_{L}$, the
Hamiltonian $H_{M}$ is in Hilbert space $\mathcal{H}_{M}$, and $H_{LM}$ gives
the coupling between the two Hilbert spaces. Note that we do not consider this
coupling to be small.

Let the dimension of $\mathcal{H}_{L}$ and $\mathcal{H}_{M}$ be $d_{L}$ and
$d_{M}$ respectively. Then, one can pick an orthonormal basis of operators
$\\{f_{i}\\}$ and $\\{g_{j}\\}$ on $\mathcal{H}_{L}$ and $\mathcal{H}_{M}$
respectively, where $1\leq i\leq d_{L}^{2}$ and $1\leq j\leq d_{M}^{2}$, and
where orthonormality is defined according to the Hilbert Schmidt inner product
given by $\braket{A,B}=\text{Tr}[A^{\dagger}B]$. One can always choose this
basis such that $f_{d_{L}^{2}}=I_{L}/\sqrt{d_{L}}$ and
$g_{d_{M}^{2}}=I_{M}/\sqrt{d_{M}}$, where $I_{M}$ and $I_{L}$ are the identity
operators on those spaces. Such a basis is required by the GKSL theorem Gorini
_et al._ (1976); Lindblad (1976); Breuer and Petruccione (2006); Rivas and
Huelga (2012). Taking the tensor product of these two basis, one can obtain an
orthonormal basis $\\{F_{k}\\}=\\{f_{i}\\}\otimes\\{g_{j}\\}$ for operators on
$\mathcal{H}_{S}$, with $F_{d_{L}^{2}d_{M}^{2}}=I_{S}/\sqrt{d}$, where
$d=d_{L}d_{M}$ is the dimension of the system Hilbert space. Without loss of
generality, the Lindblad equation [Eq. (LABEL:Lindblad_form)] written in this
basis is given by

$\frac{\partial{\rho}}{\partial
t}=i[\rho,H_{S}+H_{LS}]+\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}-1}\Gamma_{\alpha{\tilde{\alpha}}}\Big{(}F_{{\tilde{\alpha}}}\rho
F^{\dagger}_{\alpha}-\frac{\\{F^{\dagger}_{\alpha}F_{\tilde{\alpha}},\rho\\}}{2}\Big{)},$
(12)

where complete positivity of $\rho$ is preserved iff $\Gamma$ is positive
semidefinite Breuer and Petruccione (2006); Rivas and Huelga (2012). Eq. (12)
can be turned into Eq. (LABEL:Lindblad_form) by diagonalizing the matrix
$\Gamma$.

The complete positivity of RE can be checked by taking the RE to the same form
as Eq.(12) and checking if the corresponding $\Gamma$ is positive
semidefinite. To do so, let us relabel the indices so that $F_{i}=f_{i}\otimes
I_{M}/\sqrt{d_{M}}$ for $1\leq i\leq d_{L}^{2}-1$. This allows us to expand
the system operators in Eq. (8) as,

$\displaystyle S_{l}$
$\displaystyle=\sum_{\alpha=1}^{d^{2}}a_{l\alpha}F_{\alpha},\quad
S^{\dagger}_{l}=\sum_{\alpha=1}^{d^{2}}a^{\prime}_{l\alpha}F_{\alpha},$ (13)
$\displaystyle S^{(1)}_{l}$
$\displaystyle=\sum_{\alpha=1}^{d^{2}}b_{l\alpha}F_{\alpha},\quad
S^{(1)\dagger}_{l}=\sum_{\alpha=1}^{d^{2}}b^{\prime}_{l\alpha}F_{\alpha},$
$\displaystyle S^{(2)}_{l}$
$\displaystyle=\sum_{\alpha=1}^{d^{2}}c_{l\alpha}F_{\alpha},\quad
S^{(2)\dagger}_{l}=\sum_{\alpha=1}^{d^{2}}c^{\prime}_{l\alpha}F_{\alpha},$

where $a_{l\alpha}=a^{\prime}_{l\alpha}=0,\forall\quad d_{L}^{2}\leq\alpha\leq
d^{2}-1$ since $S_{l}$ and $S_{l}^{\dagger}$ are identity on
$\mathcal{H}_{M}$. Substituting Eq. (13) into Eq. (8), we obtain

$\displaystyle\frac{\partial{\rho}}{\partial
t}=i[\rho,H_{S}]-\epsilon^{2}\sum_{l}\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}\Big{\\{}a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}[F^{\dagger}_{\alpha},F_{\tilde{\alpha}}\rho]$
(14) $\displaystyle+c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}[\rho
F^{\dagger}_{\alpha},F_{\tilde{\alpha}}]+b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}[\rho
F_{\alpha}^{\dagger},F_{\tilde{\alpha}}]+a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}[F^{\dagger}_{\alpha},F_{\tilde{\alpha}}\rho]\Big{\\}}.$

Using some straightforward algebra (see Appendix A), Eq. (LABEL:RE6) can be
simplified to

$\displaystyle\frac{\partial{\rho}}{\partial t}$
$\displaystyle=i[\rho,H_{S}+H_{LS}]+\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}-1}\Gamma_{\alpha{\tilde{\alpha}}}\bigg{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}\bigg{)},$
(15)

where $\Gamma_{\alpha{\tilde{\alpha}}}$ is a $(d^{2}-1)\times(d^{2}-1)$
hermitian matrix given by

$\Gamma_{\alpha{\tilde{\alpha}}}=\epsilon^{2}\sum_{l}(a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}+c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}+b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}+a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}),$
(16)

where $x^{*}$ in Eq. (16) denotes the complex conjugate of $x$, and the
expression of $H_{LS}$ that appears in Eq. (15) can be found in Appendix A.
Since Eq. (16) is of the form in Eq. (12), the condition for preserving
complete positivity of $\rho$ is equivalent to $\Gamma$ being positive
semidefinite.

####  II.1.2 The dissipator

Let us now look at $\Gamma_{\alpha{\tilde{\alpha}}}$ from Eq. (16). Recall
from Sec. II.1.1, that if $\alpha\geq d_{L}^{2}$, then
$a_{l\alpha}=a^{\prime}_{l\alpha}=0$. Therefore,
$\Gamma_{\alpha{\tilde{\alpha}}}=0$ when $\alpha,{\tilde{\alpha}}\geq
d_{L}^{2}$, and has the following structure:

$\Gamma=\left[\begin{array}[]{ c | c }\Gamma_{\alpha,{\tilde{\alpha}}<d_{L}^{2}}&\Gamma_{\alpha<d_{L}^{2},{\tilde{\alpha}}\geq d_{L}^{2}}\\\ \hline\cr\Gamma_{\alpha\geq d_{L}^{2},{\tilde{\alpha}}<d_{L}^{2}}&0\\\ \end{array}\right].$ (17)

In general, the off-diagonal blocks
$\Gamma_{\alpha<d_{L}^{2},{\tilde{\alpha}}\geq d_{L}^{2}}$, and
$\Gamma_{\alpha\geq d_{L}^{2},{\tilde{\alpha}}<d_{L}^{2}}$ will not be
identically zero. The RE will preserve complete positivity if and only if the
matrix $\Gamma$ is positive semidefinite. To move forward we will require the
following Lemma concerning positive semidefinite matrices.

###### Lemma 1

Let $M$ be any positive semidefinite matrix, such that $M_{jj}=0$. Then,
$\forall i$ it must be the case that $M_{ij}=M_{ji}=0$. Thus, if a positive
semidefinite matrix has a zero as its $j$th diagonal element, then the entire
$j$th row and $j$th column must consist of zeros.

Proof : For a matrix $M$ to be positive semidefinite, it is necessary that
every $2\times 2$ matrix ($M^{\prime}$) of the form

$M^{\prime}=\left[\begin{array}[]{ c | c }M_{ii}&M_{ij}\\\ \hline\cr M_{ji}&M_{jj}\\\ \end{array}\right]$ (18)

is positive semidefinite. Otherwise there would exist a vector
$\ket{v^{\prime}}$ such that $\bra{v^{\prime}}M^{\prime}\ket{v^{\prime}}<0$.
Then, there would exist a vector $\ket{v}$ which is non-zero only on the $i$th
and $j$th entries, such that
$\bra{v}M\ket{v}=\bra{v^{\prime}}M^{\prime}\ket{v^{\prime}}<0$. If $M_{jj}=0$,
the eigenvalues of $M^{\prime}$ are given by solving the characteristic
equation: $\lambda^{2}-M_{ii}\lambda-M_{ij}M_{ji}=0$. Since $M$ is hermitian,
for both eigenvalues to be non-negative, it must be the case that
$M_{ij}=M_{ji}=0$.

Applying Lemma 1 to Eq. (17), we find that $\Gamma$ can be positive
semidefinite only if the off-diagonal blocks
$\Gamma_{\alpha<d_{L}^{2},{\tilde{\alpha}}\geq d_{L}^{2}}$ and
$\Gamma_{\alpha\geq d_{L}^{2},{\tilde{\alpha}}<d_{L}^{2}}$ are zero. Since
this is not generically true for RE, we have shown that the RE equation for
any generic setup where the bath only couples to a part of the system violates
complete positivity. An important exception here are situations where the RE
happens to be such that the off-diagonal blocks in $\Gamma$ are identically
zero. In such cases, the RE dissipator consists of operators of the form
$f_{i}\otimes I_{M}/\sqrt{d_{M}}$, and acts only the part of the system
connected to the bath. However, such situations can be expected to arise in
only some extremely special cases. Even in such situations, RE may or may not
preserve complete positivity, depending on
$\Gamma_{\alpha<d_{L}^{2},{\tilde{\alpha}}<d_{L}^{2}}$. Thus we have
rigorously shown that for setups where the bath couples to only a part of the
system, the RE will violate complete positivity unless its dissipator is
“local” and acts only on the part of the system that is connected to the bath.
We also see from Eq. (16), that both the positive and negative eigenvalues of
the $\Gamma$ matrix can be the same order in system-bath coupling, i.e
$O(\epsilon^{2})$.

We stress that although it is known that RE does not conserve complete
positivity, most previous works show this via specific examples (for instance,
Hartmann and Strunz (2020); Eastham  _et al._ (2016); Suárez  _et al._ (1992);
Spohn (1980)). To the best of our knowledge, we are not aware of a proof of
violation of complete positivity via the explicit, model independent
construction as shown above. In Appendix B, we give a concrete example of a
two-qubit XXZ system with the first qubit connected to the bath. We indeed
find that the matrix $\Gamma$ for this specific example has the expected
structure from Eq. (17), and negative eigenvalues.

A natural question that arises, is how one could recover complete positivity
of RE via a suitable approximation to the RE. In order to recover complete
positivity, one must change $\Gamma_{\alpha{\tilde{\alpha}}}$ to some
$\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}$ that is positive semidefinite.
By the above discussion, for any ${\tilde{\alpha}}\geq d_{L}^{2}$, this will
require either (i) making $\widetilde{\Gamma}_{\alpha\alpha}>0$,
$\forall\alpha>d_{L}^{2}$ or (ii) making
$\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}=0$ $\forall\alpha<d_{L}^{2}$. We
will see in the next section (Sec. III) that any equation with (i), will
violate local conservation laws.

##  III Lindblad descriptions obeying local conservation laws

###  III.1 Local conservation laws

Let us make precise what we mean by preservation of local conservation laws.
For our setup, since the bath only acts on $\mathcal{H}_{L}$ part of the
system, any operator on $\mathcal{H}_{M}$ commutes with the system-bath
coupling Hamiltonian $H_{SB}$ (see Eq. (10), Fig. 1). So, writing down the
Heisenberg equation of motion with respect to the full Hamiltonian $H$ [Eq.
(2)], and using Eq. (10), we see that the dynamical equation for the
expectation value of any observable $O_{M}$ on $\mathcal{H}_{M}$ is given by

$\frac{d}{dt}\braket{I_{L}\otimes O_{M}}=-i\braket{[I_{L}\otimes
O_{M},H_{S}]}$ (19)

where $\braket{X}=\text{Tr}[X\rho]$. Any effective QME obtained by integrating
out the bath should satisfy this property. We call QMEs satisfying this
property as ones preserving local conservation laws. The justification for
this name becomes clear if we look at an operator in $\mathcal{H}_{M}$ that
would remain conserved if there is coupling with $\mathcal{H}_{L}$, i.e, if
$H_{LM}$ in Eq.(11) is zero. One such operator is the Hamiltonian $H_{M}$. The
dynamical equation for expectation value of $H_{M}$ gives the energy
continuity equation

$\displaystyle\frac{d}{dt}\braket{H_{M}}=J_{L\to M},~{}~{}J_{L\to
M}=-i\braket{[H_{M},H_{LM}]}.$ (20)

Here, $J_{L\to M}$ can be interpreted as the energy current from the region
$L$ to the region $M$ (see Fig. 1). In steady state, the rate of change of any
system operator is zero. From above equation, this gives, $J_{L\to M}=0$ in
steady state. Thus, steady state energy current inside the system is zero.
This is a statement of local conservation of energy and is one of the
fundamental physical requirements for a system coupled to a single bath that
follows from the more general requirement Eq. (19).

Importantly, the RE, i.e, Eq. (5), can be shown to satisfy Eq. (19) Tupkary
_et al._ (2022) and thereby preserves local conservation laws. We can write
any QME to leading order in system bath coupling as

$\frac{\partial{\rho}}{\partial
t}=\mathcal{L}_{0}(\rho)+\epsilon^{2}\mathcal{L}_{2}(\rho),$ (21)

where $\mathcal{L}_{0}(\rho)=i[\rho,H_{S}]$ and $\mathcal{L}_{2}(\rho)$
contains both the dissipator and the Lamb-shift Hamiltonian. Computing the
left hand size of Eq. (19) using Eq. (21), and comparing with the right hand
size of Eq. (19), we obtain Tupkary  _et al._ (2022),

$\text{Tr}[(I_{M}\otimes O_{M})\mathcal{L}_{2}(\rho)]=0.$ (22)

This is a necessary condition for satisfying local conservation laws. If we
now further restrict the QME to be of Lindblad form, i.e, of the form Eq.
(LABEL:Lindblad_form), thereby respecting complete positivity, we obtain the
following theorem.

###### Theorem 2

Any QME of Lindblad form Eq. (LABEL:Lindblad_form) (thereby satisfying
complete positivity) that also satisfies local conservation laws must have the
Lindblad operators and the Lamb-shift Hamiltonian acting only on the part of
the system connected to the bath. That is, $L_{\lambda}=L_{\lambda}^{L}\otimes
I_{M}$, $H_{LS}=H_{LS}^{L}\otimes I_{M}$, where $L_{\lambda}^{L}$,
$H_{LS}^{L}$ act only on the Hilbert space $\mathcal{H}_{L}$ which is coupled
to the bath by system-bath coupling Hamiltonian $H_{SB}$, and $I_{M}$ is the
identity on the remaining of the system Hilbert space $\mathcal{H}_{M}$.

In the next subsection, we give the proof of this theorem.

###  III.2 Proof of theorem 2

We start by writing the most general Lindblad equation in the basis of Sec.
II.1.1,

$\frac{\partial{\rho}}{\partial
t}=i[\rho,H_{S}+\widetilde{H}_{LS}]+\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}-1}\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}\Big{(}F_{{\tilde{\alpha}}}\rho
F^{\dagger}_{\alpha}-\frac{\\{F^{\dagger}_{\alpha}F_{\tilde{\alpha}},\rho\\}}{2}\Big{)}$
(23)

where $\widetilde{H}_{LS}$ is some Lamb shift Hamiltonian, and
$\widetilde{\Gamma}$ is a positive semidefinite matrix. As mentioned before,
this form can be reduced to the form of Eq.(LABEL:Lindblad_form) by
transforming to a basis where matrix $\widetilde{\Gamma}$ is diagonal. So, it
suffices to work with this form. Since $F_{i}=f_{i}\otimes I_{M}$ for $1\leq
i<d_{L}^{2}$, the condition for the Lindblad operators to act only on
$\mathcal{H}_{L}$ then translates to the matrix $\widetilde{\Gamma}$ being of
the form

$\widetilde{\Gamma}=\left[\begin{array}[]{ c | c }\widetilde{\Gamma}_{\alpha,{\tilde{\alpha}}<d_{L}^{2}}&0\\\ \hline\cr 0&0\\\ \end{array}\right].$ (24)

In the following we prove that in order to preserve local conservation laws
the matrix $\widetilde{\Gamma}$ must be of this form.

####  III.2.1 The restriction on $\widetilde{\Gamma}$

Writing down the evolution of any observable $I_{L}\otimes O_{M}$ for Eq.
(23), we have

$\displaystyle\frac{d\braket{I_{L}\otimes O_{M}}}{dt}=-i\braket{[I_{L}\otimes
O_{M},H_{S}]}-i\braket{[I_{L}\otimes O_{M},\widetilde{H}_{LS}]}$ (25)
$\displaystyle+\sum_{\alpha=1}^{d_{L}^{2}-1}\sum_{{\tilde{\alpha}}=d_{L}^{2}}^{d^{2}-1}\frac{\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}}{2}\braket{[I_{L}\otimes
O_{M},F^{\dagger}_{\alpha}F_{\tilde{\alpha}}]}$
$\displaystyle+\sum_{\alpha=d_{L}^{2}}^{d^{2}-1}\sum_{{\tilde{\alpha}}=1}^{d_{L}^{2}-1}\frac{\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}}{2}\braket{[F^{\dagger}_{\alpha}F_{\tilde{\alpha}},I_{L}\otimes
O_{M}]}$
$\displaystyle+\sum_{\alpha,{\tilde{\alpha}}=d_{L}^{2}}^{d^{2}-1}\Big{(}\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}\Big{(}\braket{F^{\dagger}_{\alpha}(I_{L}\otimes
O_{M})F_{\tilde{\alpha}}}-\frac{1}{2}\braket{(I_{L}\otimes
O_{M})F_{\alpha}^{\dagger}F_{\tilde{\alpha}}}$
$\displaystyle-\frac{1}{2}\braket{F_{\alpha}^{\dagger}F_{\tilde{\alpha}}(I_{L}\otimes
O_{M})}\Big{)}$

where $\braket{X}$ denotes the expectation value of $X$ given by
$\text{Tr}[X\rho]$, and we use the fact that $[F_{i},I_{L}\otimes O_{M}]=0$
when $i<d_{L}^{2}$. We can combine all commutators in above equation into a
single commutator with an effective Lamb-shift like Hamiltonian, which we
denote $\widetilde{H}^{(2)}_{LS}$. This gives

$\displaystyle\frac{d\braket{I_{L}\otimes O_{M}}}{dt}=-i\braket{[I_{L}\otimes
O_{M},H_{S}+\widetilde{H}^{(2)}_{LS}]}$ (26)
$\displaystyle+\sum_{\alpha,{\tilde{\alpha}}=d_{L}^{2}}^{d^{2}-1}\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}\Big{(}\braket{F^{\dagger}_{\alpha}(I_{L}\otimes
O_{M})F_{\tilde{\alpha}}}-\frac{1}{2}\braket{(I_{L}\otimes
O_{M})F_{\alpha}^{\dagger}F_{\tilde{\alpha}}}$
$\displaystyle-\frac{1}{2}\braket{F_{\alpha}^{\dagger}F_{\tilde{\alpha}}(I_{L}\otimes
O_{M})}\Big{)}$

where $\widetilde{H}^{(2)}_{LS}$ is some hermitian operator. Comparing with
Eq. (19), and noting that $\operatorname{Tr}[M\rho]=0\quad\forall\rho$ implies
$M=0$, we obtain the condition for satisfying Eq. (19) as

$\displaystyle-i[I_{L}\otimes
O_{M},\widetilde{H}^{(2)}_{LS}]+\sum_{\alpha,{\tilde{\alpha}}=d_{L}^{2}}^{d^{2}-1}\widetilde{\Gamma}_{\alpha{\tilde{\alpha}}}\Big{(}F^{\dagger}_{\alpha}(I_{L}\otimes
O_{M})F_{\tilde{\alpha}}$ (27) $\displaystyle-\frac{1}{2}(I_{L}\otimes
O_{M})F_{\alpha}^{\dagger}F_{\tilde{\alpha}}-\frac{1}{2}F_{\alpha}^{\dagger}F_{\tilde{\alpha}}(I_{L}\otimes
O_{M})\Big{)}=0$ $\displaystyle\hskip 190.00029pt\forall O_{M}.$

We will now relabel the indices $\alpha,{\tilde{\alpha}}$ for the sake of
convenience. Recall from Sec. II.1.1, that $\alpha,{\tilde{\alpha}}$ which
appear in the above expression can be equivalently written as
$\alpha\rightarrow(\alpha_{L},\alpha_{M}),{\tilde{\alpha}}\rightarrow({\tilde{\alpha}}_{L},{\tilde{\alpha}}_{M})$,
and vice-versa, where $F_{\alpha}=f_{\alpha_{L}}\otimes g_{\alpha_{M}}$.
Therefore, $\sum_{\alpha=1}^{d^{2}-1}$ is equivalent to
$\sum_{\alpha_{M}=1}^{d^{2}_{M}-1}\sum_{\alpha_{L}=1}^{d_{L}^{2}}$. We can
also expand $\widetilde{H}^{(2)}_{LS}$ as

$\widetilde{H}^{(2)}_{LS}=\sum_{\alpha=1}^{d^{2}}\nu_{\alpha}F_{\alpha}=\sum_{\alpha_{L}=1}^{d^{2}_{L}}\sum_{\alpha_{M}=1}^{d^{2}_{M}}\nu_{\alpha_{L},\alpha_{M}}(f_{\alpha_{L}}\otimes
g_{\alpha_{M}}).$ (28)

This basis aids in taking a partial trace over the $L$ part of the system.
Performing this partial trace and using the orthonormality of $f_{i}$
operators, we can rewrite Eq. (27) as

$\displaystyle-i\sum_{\alpha_{M}=1}^{d_{M}^{2}}\nu_{d_{L}^{2},\alpha_{M}}[O_{M},g_{\alpha_{M}}]+\sum_{\alpha_{M},{\tilde{\alpha}}_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}$
(29)
$\displaystyle\Big{(}g_{\alpha_{M}}^{\dagger}O_{M}g_{{\tilde{\alpha}}_{M}}-\frac{1}{2}O_{M}g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}-\frac{1}{2}g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}O_{M}\Big{)}=0,$
$\displaystyle\hskip 190.00029pt\forall O_{M}.$

where

$\displaystyle\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}$
$\displaystyle=\sum_{\alpha_{L},{\tilde{\alpha}}_{L}=1}^{d_{L}^{2}}\delta_{\alpha_{L},{\tilde{\alpha}}_{L}}\widetilde{\Gamma}_{(\alpha_{L},\alpha_{M}),({\tilde{\alpha}}_{L},{\tilde{\alpha}}_{M})}$
(30)

We show in Appendix C that Eq. (LABEL:cont_conjecture_2) implies

$\sum_{\alpha_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},\alpha_{M}}=\sum_{\alpha_{M}=1}^{d_{M}^{2}-1}\sum_{\alpha_{L}=1}^{d_{L}^{2}}\widetilde{\Gamma}_{(\alpha_{L},\alpha_{M}),(\alpha_{L},\alpha_{M})}=0$
(31)

Since we require $\widetilde{\Gamma}$ to be positive semidefinite, it cannot
have negative values on the diagonal. Therefore, Eq. (31) implies
$\widetilde{\Gamma}_{(\alpha_{L},\alpha_{M}),(\alpha_{L},\alpha_{M})}=0$ for
$1\leq\alpha_{L}\leq d_{L}^{2}$ and $1\leq\alpha_{M}\leq d_{M}^{2}-1$ .
Equivalently, $\widetilde{\Gamma}_{\alpha\alpha}=0$ for $\alpha\geq
d_{L}^{2}$. Applying Lemma 1 to this case, we see that
$\widetilde{\Gamma}_{\alpha,{\tilde{\alpha}}}$ can be non-zero only when both
$\alpha,{\tilde{\alpha}}<d_{L}^{2}$, and therefore $\widetilde{\Gamma}$ is a
matrix of the form in Eq. (24). This concludes the first part of the proof.

####  III.2.2 The restriction on $\widetilde{H}_{LS}$

Given this structure for $\widetilde{\Gamma}$ from Eq. (24), we now
investigate the restrictions on the Lamb shift Hamiltonian
$\widetilde{H}_{LS}$. Since $\widetilde{\Gamma}$ has to obey Eq. (24), we find
that
$\widetilde{H}_{LS}=\widetilde{H}_{LS}^{(2)}=\sum_{\alpha_{L},\alpha_{M}}\nu_{\alpha_{L},\alpha_{M}}(f_{\alpha_{L}}\otimes
g_{\alpha_{M}})$ in Eq. (27). Then, our condition for satisfying local
conservation laws from Eq. (27) is given by,

$-i\sum_{\alpha_{L}=1}^{d^{2}_{L}}\sum_{\alpha_{M}=1}^{d^{2}_{M}}\nu_{\alpha_{L},\alpha_{M}}[I_{L}\otimes
O_{M},(f_{\alpha_{L}}\otimes g_{\alpha_{M}})]=0\quad\forall O_{M}.$ (32)

which further implies

$-i\sum_{\alpha_{L}=1}^{d^{2}_{L}}\sum_{\alpha_{M}=1}^{d^{2}_{M}}\nu_{\alpha_{L},\alpha_{M}}f_{\alpha_{L}}\otimes[O_{M},g_{\alpha_{M}}]=0\quad\forall
O_{M}.$ (33)

Multiplying both sides by $f^{\dagger}_{\alpha_{L}}\otimes I_{M}$, and tracing
out the $L$ part of the system, we obtain

$-i\sum_{\alpha_{M}=1}^{d^{2}_{M}}\nu_{\alpha_{L},\alpha_{M}}[O_{M},g_{\alpha_{M}}]=0\quad\forall\alpha_{L},O_{M}.$
(34)

This can happen only if

$\sum_{\alpha_{M}=1}^{d^{2}_{M}}\nu_{\alpha_{L},\alpha_{M}}g_{\alpha_{M}}\propto
I_{M}\quad\forall\alpha_{L}$ (35)

Using Eq. (35) in Eq. (28), and recalling that
$\widetilde{H}_{LS}=\widetilde{H}^{(2)}_{LS}$, we obtain
$\widetilde{H}_{LS}=\widetilde{H}^{(L)}_{LS}\otimes I_{M}$. This concludes the
second part of the proof.

###  III.3 Remarks on Theorem 2

Theorem 2 says that for a QME to preserve complete positivity and obey local
conservation laws, both the Lamb shift Hamitonian and the dissipator must only
act on the part of the system connected to the bath. Such a Lindblad equation
is often termed a ‘local Lindblad equation’. Theorem 2 thus says that only
local Lindblad equations are consistent with local conservation laws. The RE
preserves local conservation laws without having a local dissipator, but it
does so at the cost of losing complete positivity. Any global form of Lindblad
equation, like the eigenbasis Lindblad equation Carmichael (2002); Breuer and
Petruccione (2006) and the universal Lindblad Equation Nathan and Rudner
(2020), violates local conservation laws, while preserving complete
positivity. One main reason such global forms of Lindblad equations are often
derived under various approximations is that they can be proven to show
thermalization. However, general statements related to thermalization in local
Lindblad equations are usually difficult to make. In the next section, we
discuss a numerical technique which allows to study whether, in a given set-
up, a QME consistent with Theorem 2 is possible such that it also shows
thermalization.

##  IV On the possibility of thermalization with local dissipators

###  IV.1 Condition for satisfying thermalization

We start by making precise what we mean by thermalization. Going back to the
form of the QME in Eq. (21), the setup is said to show thermalization if

$\displaystyle\lim_{\epsilon\to
0}\left(\lim_{t\to\infty}e^{t(\mathcal{L}_{0}+\epsilon^{2}\mathcal{L}_{2})}\rho(0)\right)=\rho_{th},~{}~{}\rho_{th}=\frac{e^{-\beta
H_{S}}}{\text{Tr}[e^{-\beta H_{S}}]},$ (36)

irrespective of the initial state of the system $\rho(0)$. Physically, it
means that if the system is weakly coupled to a thermal bath for a long time,
and then the system-bath coupling is slowly switched-off, the state of the
system will be the Gibbs state at the temperature of the bath, irrespective of
the system’s initial state. If there is no explicit time-dependence in the
Hamiltonian, as in Eq. (2), this statement can be proven starting with the
initial state of the full set-up in Eq. (3), and assuming that the steady
state is unique Tupkary  _et al._ (2022). Given a QME of the form in Eq.(21),
in Ref. Tupkary  _et al._ (2022), we showed that the following condition needs
to be satisfied in order to guarantee thermalization,

$\bra{E_{i}}\mathcal{L}_{2}(\rho_{th})\ket{E_{i}}=0\quad\forall~{}i.$ (37)

where $\ket{E_{i}}$ is the eigenvector of the system Hamiltonian $H_{S}$ with
eigenvalue $E_{i}$. The derivation of this condition can also be found in
Appendix D. Given a system Hamiltonian $H_{S}$, the inverse temperature
$\beta$ of the bath and a partition of the system Hilbert space into the part
$\mathcal{H}_{L}$ which is attached to a bath and the remainder of Hilbert
space $H_{M}$ (see Fig.1), we would like to find a QME respecting the
restrictions in theorem 2 and satisfying Eq.(37). As we will see later by
example, such a QME is not guaranteed to be possible. In the next subsection,
we provide a numerical way to conclusively check if such a QME is possible in
a given setup.

###  IV.2 The thermalization optimization problem

For the setup shown in Fig. 1, the most general form of QME respecting the
restrictions in theorem 2 for satisfying complete positivity and local
conservation laws is

$\displaystyle\frac{\partial{\rho}}{\partial
t}=i[\rho,H_{S}+\epsilon^{2}H^{(L)}_{LS}\otimes
I_{M}]+\epsilon^{2}\sum_{\alpha_{L},{\tilde{\alpha}}_{L}=1}^{d_{L}^{2}-1}\Gamma^{(L)}_{\alpha_{L},{\tilde{\alpha}}_{L}}$
(38) $\displaystyle\bigg{(}(f_{\tilde{\alpha}_{L}}\otimes
I_{M})\rho(f_{\alpha_{L}}\otimes I_{M})^{\dagger}$
$\displaystyle-\frac{\\{(f_{\alpha_{L}}\otimes
I_{M})^{\dagger}(f_{\tilde{\alpha}_{L}}\otimes I_{M}),\rho\\}}{2}\bigg{)}.$

Here $H^{(L)}_{LS}$ is a Lamb shift Hamiltonian that acts on the
$\mathcal{H}_{L}$ part of the system, $\Gamma^{(L)}$ is
$(d_{L}^{2}-1)\otimes(d_{L}^{2}-1)$ matrix that must be positive semidefinite,
$d_{L}$ is the dimension of the Hilbert space $\mathcal{H}_{L}$ that is
directly coupled with the bath. We include the factor of $\epsilon^{2}$ in
front of $\Gamma^{(L)}$ and $H^{(L)}_{LS}$ explicitly. The system Hamiltonian
$H_{S}$ and the Hilbert space $\mathcal{H}_{L}$ are assumed to be given. The
task is then to find $H^{(L)}_{LS}$ and $\Gamma^{(L)}$, such that Eq.(37) is
satisfied to a given precision, for a given inverse temperature $\beta$. To
this end, we introduce the quantity

$\tau=\sum_{i}\left|\bra{E_{i}}\mathcal{L}_{2}(\rho_{th})\ket{E_{i}}\right|,$
(39)

where $\mathcal{L}_{2}$ consist of all terms in Eq.
(LABEL:eq:local_lindblad_for_TOP) which are multiplied by $\epsilon^{2}$, i.e
all terms except the commutator with the system Hamiltonian. Then, we can cast
the task in terms of an optimization problem given by :

minimize : $\displaystyle\tau~{}\textrm{by varying
}H^{(L)}_{LS},~{}\Gamma^{(L)},$ (40) subject to : $\displaystyle
H^{(L)}_{LS}\text{ is hermitian,}$
$\displaystyle\text{Tr}(\Gamma^{(L)})=1,~{}~{}\Gamma^{(L)}\geq 0,$

where we use $\Gamma^{(L)}\geq 0$ to denote $\Gamma^{(L)}$ being positive
semidefinite. The condition $\text{Tr}(\Gamma^{(L)})=1$ is imposed to avoid
the trivial solution $H^{(L)}_{LS},\Gamma^{(L)}=0$, which trivially gives the
global minimum $\tau=0$. Since we want $\Gamma^{(L)}$ to be a non-zero
positive semi-definite matrix, it must have a positive trace. We have fixed
that trace to one in some arbitrarily chosen energy unit. This does not cause
loss of generality since the strength of system-bath coupling is explicitly
governed by $\epsilon^{2}$ in Eq.(LABEL:eq:local_lindblad_for_TOP). We call
the optimization problem in Eq. (40) the “thermalization optimization
problem”(TOP). Let $\tau_{opt}$ be the optimal value obtained from solving
TOP. Given a tolerance $\delta$,

$\displaystyle\begin{array}[]{cc}\textrm{if }\tau_{opt}<\delta,&\textrm{ the
desired QME is possible},\\\ \textrm{else, }&\textrm{ the desired QME is
impossible},\end{array}$ (43)

up to the precision $\delta$.

Most interestingly, we find that the TOP in Eq. (40) can be written as a SDP.
The background and theoretical framework of SDP is discussed in Appendix E.1.
The TOP is proven to be a SDP in Appendix E.2. In particular, any choice of
$H^{(L)}_{LS}$ and $\Gamma^{(L)}$ in Eq. (40) can be used to obtain an upper
bound on $\tau_{opt}$. Then, the theoretical framework of SDP can be used to
construct a “dual” problem to the optimization problem in Eq. (40). This dual
problem can then be used to obtain a lower bound on $\tau_{opt}$. If the lower
bound and upper bound match, then this guarantees that one has found the
global optimal value of $\tau$ . Our above described approach is transparent,
and simple to use, since the optimization problem in Eq. (40) can be directly
put into standard packages for disciplined convex optimization like the CVX
MATLAB package Grant and Boyd (2014). In particular, CVX itself automatically
constructs the dual problem, outputs $\tau_{opt}$, and gives one choice of
$H^{(L)}_{LS}$, and $\Gamma^{(L)}$ which yields the output value of
$\tau_{opt}$. Thus, if $\tau_{opt}<\delta$, it not only says that the desired
type of QME is possible but also it outputs one possible candidate for such a
QME. If $\tau_{opt}\geq\delta$, the desired type of QME is impossible.

We would like to point out here that, in a microscopic derivation, given the
temperature of the baths, $H^{(L)}_{LS}$ and $\Gamma^{(L)}$ would depend only
on the bath spectral functions and system-bath coupling Hamiltonian. So, the
TOP can be thought of as varying over all possible bath spectral functions and
system-bath coupling Hamiltonians to find $\tau_{opt}$. Thus, if
$\tau_{opt}\geq\delta$, we can conclusively say that, for the chosen system
parameters and temperature, under no choice of bath spectral function and
system-bath coupling Hamiltonian, can a Markovian QME be derived which
simultaneously satisfies complete positivity, local conservation laws and
shows thermalization up to the chosen precision. In the next subsection, we
look at the TOP in an open XXZ qubit chain.

###  IV.3 Open XXZ qubit chain as an example

We study the possibility of having a Lindblad description satisfying local
conservation laws and showing thermalization in an open XXZ qubit chain system
with some of the qubits attached to baths. The system Hamiltonian for this
setup is given by

$\displaystyle
H_{S}=\sum_{\ell=1}^{N}\frac{\omega^{(\ell)}_{0}}{2}{\sigma}_{z}^{(\ell)}-\sum_{\ell=1}^{N-1}g_{\ell}\big{(}$
$\displaystyle{\sigma}_{x}^{(\ell)}{\sigma}_{x}^{(\ell+1)}+{\sigma}_{y}^{(\ell)}{\sigma}_{y}^{(\ell+1)}$
(44)
$\displaystyle+\Delta_{\ell}{\sigma}_{z}^{(\ell)}{\sigma}_{z}^{(\ell+1)}\big{)},$

where $\sigma^{(\ell)}_{x,y,z}$ denotes the Pauli matrices acting on the
$\ell^{\text{th}}$ qubit, and $\omega^{(\ell)}_{0}$, $g_{\ell}$, and
$g_{\ell}\Delta_{\ell}$ represent the magnetic field, the overall qubit-qubit
coupling strength and the anisotropy respectively. The first $N_{L}$ qubits
are attached to a bath, while the remaining $N_{M}=N-N_{L}$ qubits are not
attached to any bath. We use the formalism of Sec. IV.2 to investigate
thermalization in this set-up for various values of $N_{L}$ and $N_{M}$. In
order to do so, we first need to construct the basis for operators in Hilbert
space of the first $N_{L}$ qubits, so that the most general form of the
desired QME can be written as in Eq. (LABEL:eq:local_lindblad_for_TOP). For
the $\ell$th qubit, we choose the basis
$\\{-\sigma^{(\ell)}_{z}/\sqrt{2},\sigma^{(\ell)}_{+},\sigma^{(\ell)}_{-},I^{(\ell)}_{2}/\sqrt{2}\\}$,
where $\sigma^{(\ell)}_{+}=(\sigma^{(\ell)}_{x}+i\sigma^{(\ell)}_{y})/2$,
$\sigma^{(\ell)}_{-}=(\sigma^{(\ell)}_{x}-i\sigma^{(\ell)}_{y})/2$, and
$I^{(\ell)}_{2}$ is the identity operator for the qubit Hilbert space. The
basis for the first $N_{L}$ qubits is obtained by direct product of the basis
of each of the qubits. We construct the TOP [Eq.(40)] in this basis, which we
then directly input in the CVX MATLAB package to obtain $\tau_{opt}$. We set
an ad-hoc value of the tolerance $\delta=10^{-6}$ [see Eq.(43)]. In a typical
calculation from a weak system-bath coupling QME, usually the error due to
neglecting higher order terms would be larger than such a low value of
tolerance.

Figure 2: $\tau_{opt}$ vs $g$, for $N_{L}=1$, with $\omega^{(\ell)}_{0}=1$,
$\Delta_{\ell}=1$, $\beta=1$ and $g_{\ell}=g$ for all $\ell$. The tolerance
chosen is $\delta=10^{-6}$. We find that $\tau_{opt}\gg\delta$, conclusively
showing that, for such setups, no QME can simultaneously preserve complete
positivity, obey local conservation laws, and show thermalization up to the
precision set by the tolerance.  Figure 3: $\tau_{opt}$ vs $\beta$, for
$N_{L}=1$, with $\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$, $g_{\ell}=0.1$ for
all $\ell$. The tolerance chosen is $\delta=10^{-6}$, and is plotted as the
dashed horizontal line. We find that $\tau_{opt}\gg\delta$, indicating that,
for such setups, it is not possible to have a QME simultaneously preserving
complete positivity, obeying local conservation laws, and showing
thermalization up to the precision set by the tolerance.

####  IV.3.1 Single qubit attached to bath

First we consider the case where the first qubit of the chain is coupled to a
bath, so, $N_{L}=1$. In Fig. (2), we plot $\tau_{opt}$ vs $g$ for
$N_{M}=1,2,3$, fixing $\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$ for all
$\ell$ and $\beta=1$. We find that $\tau_{opt}\gg\delta$ when $N_{L}=1$, and
$N_{M}=1,2,3$. Thus, no Markovian QME can simultaneously preserve complete
positivity, obey local conservation laws, and satisfy thermalization for such
setups in the chosen range of parameters. This explicit example also directly
rules out the possibility of having a general form of Markovian QME that is
guaranteed to meet all the fundamental requirements. We see from Fig. (2) that
$\tau_{opt}$ increases with $g$. This is consistent with previous results
showing that such local Lindblad equations become a good description when the
coupling between the system qubits are weak. Additionally, we see that for a
given value of $g$, $\tau_{opt}$ decreases as $N_{M}$ increases. This,
interestingly, seems to indicate that for long chains, local Lindblad with
dissipator acting on only one qubit might be able to describe the thermal
state up to a reasonably good precision. However, more detailed studies are
required to make any conclusive statement in this direction.

Intuitively, memory effects, and hence Markovianity of open system dynamics,
depend on temperature. So, we might expect that at a different temperature,
$\tau_{opt}$ might be smaller than $\delta$. In Fig. (3), we plot $\tau_{opt}$
vs $\beta$ for $N_{L}=1$, fixing $\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$,
$g_{\ell}=0.1$ for all $\ell$. We find that $\tau_{opt}\gg\delta$, for almost
the entire range of $\beta$ chosen, showing that for these parameters, no
Markovian QME can simultaneously preserve complete positivity, obey local
conservation laws, and satisfy thermalization. But, we see some interesting
features. Firstly, as before, we see that $\tau_{opt}$ decreases as $N_{M}$
increases. Secondly, we see that $\tau_{opt}$ varies non-monotonically with
$\beta$. At very low temperatures, $\tau_{opt}$ decreases tends to decay below
$\delta$. In Fig. (3), for $N_{M}=4$ and $\beta=10$, $\tau_{opt}<\delta$. At
high temperatures also, $\tau_{opt}$ decreases. This suggests that at very low
and very high temperatures, it is possible to obtain a local Lindblad equation
that shows thermalization. That this is indeed true can be checked
independently. At such extremes of temperatures, for at least some choices of
baths and system-bath couplings, local Linblad equations can be
microscopically derived Walls and Milburn (2008); Becker  _et al._ (2021).
This explains the non-monotonic dependence of $\tau_{opt}$ on $\beta$. Since
very high and very low temperatures can allow for a local Lindblad
description, it is then intuitive that departure from such behavior is maximum
when $\beta$ is of the order of system time scales. Indeed it is close to such
values of $\beta$, i.e, $\beta\sim\omega_{0}^{(\ell)}=1$, that the highest
value of $\tau_{opt}$ is seen in Fig. (3).

We would like to point out here that, $\tau_{opt}\geq\delta$ does not mean the
system coupled to a thermal bath is unable to thermalize for such parameters.
Unless in extremely special cases, it is almost always possible to find a bath
spectral function and system-bath coupling which ensure thermalization. So
$\tau_{opt}\geq\delta$ instead means that, for those parameters, the dynamics
of approach to thermal state cannot be governed by a completely positive
Markovian QME preserving local conservation laws. The dynamics then must have
some non-Markovian character. In fact, as we see from above example,
$\tau_{opt}$ gives an estimate of how close to Markovianity the open system
dynamics can be for the chosen parameters.

Next, we discuss the case where first two qubits are attached to the bath and
highlight the drastic difference observed for similar choice of parameters.

Figure 4: $\tau_{opt}$ vs $g$, for $N_{L}=2$, with $\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$, $g_{\ell}=g$ for all $\ell$ and $\beta=1$. The tolerance chosen is $\delta=10^{-6}$, and is plotted as the dashed horizontal line. We find that $\tau_{opt}\ll\delta$ for smaller values of $g$, indicating that, for such setups, it is possible to have a QME simultaneously preserving complete positivity, obeying local conservation laws, and showing thermalization up to the precision set by the tolerance.  Figure 5: $\tau_{opt}$ vs $\beta$, for $N_{L}=2$, with $\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$, $g_{\ell}=0.1$ for all $\ell$. The tolerance chosen is $\delta=10^{-6}$, and is plotted as the dashed horizontal line. We find that $\tau_{opt}\ll\delta$ for smaller values of $g$, indicating that, for such setups, it is possible to have a QME simultaneously preserving complete positivity, obeying local conservation laws, and showing thermalization up to the precision set by the tolerance. The values in the figure less than $10^{-12}$ are below the numerical precision of CVX Matlab package.  0.00033 | 0 | 0 | -0.00123 | 0 | 0 | 0.00033 | 0 | 0 | 0.00019 | 0 | 0 | 0.00001 | 0 | 0
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
0 | 0.00196 | 0 | 0 | -0.00179 | 0 | 0 | -0.03638 | 0 | 0 | 0 | 0 | 0 | -0.00002 | 0
0 | 0 | 0.00081 | 0 | 0 | 0 | 0 | 0 | -0.00081 | 0 | 0 | 0.01307 | 0 | 0 | -0.00021
-0.00123 | 0 | 0 | 0.02611 | 0 | 0 | -0.00282 | 0 | 0 | 0.00157 | 0 | 0 | -0.00143 | 0 | 0
0 | -0.00179 | 0 | 0 | 0.00172 | 0 | 0 | 0.03491 | 0 | 0 | 0 | 0 | 0 | 0.00011 | 0
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
0.00033 | 0 | 0 | -0.00282 | 0 | 0 | 0.00045 | 0 | 0 | 0.00002 | 0 | 0 | 0.00011 | 0 | 0
0 | -0.03638 | 0 | 0 | 0.03491 | 0 | 0 | 0.70685 | 0 | 0 | 0 | 0 | 0 | 0.00209 | 0
0 | 0 | -0.00081 | 0 | 0 | 0 | 0 | 0 | 0.00082 | 0 | 0 | -0.01311 | 0 | 0 | 0.00021
0.00019 | 0 | 0 | 0.00157 | 0 | 0 | 0.00002 | 0 | 0 | 0.00035 | 0 | 0 | -0.00014 | 0 | 0
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
0 | 0 | 0.01307 | 0 | 0 | 0 | 0 | 0 | -0.01311 | 0 | 0 | 0.26032 | 0 | 0 | -0.00208
0.00001 | 0 | 0 | -0.00143 | 0 | 0 | 0.00011 | 0 | 0 | -0.00014 | 0 | 0 | 0.00009 | 0 | 0
0 | -0.00002 | 0 | 0 | 0.00011 | 0 | 0 | 0.00209 | 0 | 0 | 0 | 0 | 0 | 0.0001 | 0
0 | 0 | -0.00021 | 0 | 0 | 0 | 0 | 0 | 0.00021 | 0 | 0 | -0.00208 | 0 | 0 | 0.00009
Table 1:  The $\Gamma^{(L)}$ obtained from CVX for $N_{L}=2$, $N_{M}=4$,
$\omega_{0}^{(\ell)}=1$, $g_{\ell}=0.1$, $\Delta_{\ell}=1$, $\beta=1$, i.e,
the values used to compute Fig. 6. Every entry is rounded to 5 digits after
the decimal point for convenience of representation. The corresponding
$H^{(L)}_{LS}$ is given to be zero. This $\Gamma^{(L)}$ satisfies complete
positivity, local conservation laws and thermalization up to a precision of
$\delta=10^{-6}$ for the given choice of parameters.

####  IV.3.2 Two qubits attached to bath

In Fig. (4), we plot $\tau_{opt}$ vs $g_{\ell}=g$ for $N_{L}=2$, i.e, first
two qubits attached to the bath, and $N_{M}=1,2,3,4$, taking fixed values of
$\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$ and $\beta=1$. These parameters are
the same as in Fig. 2 and $g$ is varied over the same range. Quite remarkably,
in stark contrast to Fig. 2, we find that in this case $\tau_{opt}\ll\delta$
over a considerable range for $g<1$. We then look at the behavior of
$\tau_{opt}$ versus $\beta$, fixing $\omega^{(\ell)}_{0}=1$,
$\Delta_{\ell}=1$, $g_{\ell}=g=0.1$. This is shown in Fig. 5. For
$N_{M}>N_{L}$, we again see the non-monotonic behavior. However, in stark
contrast to Fig. 3, we find that over the entire chosen range of $\beta$
$\tau_{opt}\ll\delta$. Thus, for $N_{L}=2$, over a considerable range of
parameters, a QME that simultaneously preserves complete positivity, obeys
local conservation laws, and satisfies thermalization up to the given
precision is possible. This is a highly non-trivial result.

Previously, local two-qubit Lindblad dissipators have been used to study
energy transport in XXZ-type qubit chains (for example, Refs. Prosen and
Žnidarič (2009); Mendoza-Arenas  _et al._ (2015)). However, those local two-
qubit Lindblad operators were constructed so as to thermalize the two qubits
only, in absence of coupling to the rest of the chain. Such Lindblad
description is not guaranteed to thermalize the whole chain to the given
inverse temperature $\beta$ of the bath Žnidarič  _et al._ (2010). Our result
here shows that it is possible to have a two-qubit local Lindblad description
that can thermalize the full chain to the given temperature of the bath to a
good approximation.

As mentioned before, CVX also outputs a possible choice of $\Gamma^{(L)}$ and
$H^{(L)}_{LS}$ matrices corresponding to $\tau_{opt}$. So, when
$\tau_{opt}<\delta$, we get one possible candidate for the desired type of
QME. For our choice of parameters, we find that CVX always outputs
$H^{(L)}_{LS}=0$, and a non-trivial value of $\Gamma^{(L)}$ that would be hard
to guess otherwise. In Table. 1, we demonstrate the $\Gamma^{(L)}$ obtained
for $N_{L}=2$, $N_{M}=4$, $\omega_{0}^{(\ell)}=1$, $g_{\ell}=0.1$,
$\Delta_{\ell}=1$, $\beta=1$. The $\Gamma^{(L)}$ matrix corresponds to the
basis of operators $\\{F_{k}\\}$ chosen as

$F_{k}=f_{\lceil k/4\rceil}\otimes f_{k(\text{mod }4)}\otimes I_{M}$ (45)

where
$\\{f_{i}\\}=\\{-\sigma_{z}/\sqrt{2},\sigma_{-},\sigma_{+},I_{2}/\sqrt{2}\\}$,
$\lceil k/4\rceil$ denotes the nearest integer greater than or equal to $k/4$,
and $k(\text{mod }4)$ denotes the value of $k$ modulo $4$, and $k$ goes from
$1$ to $15$. We also note that the exact values of $\Gamma^{(L)}$ and
$H^{(L)}_{LS}$ computed by CVX may depend on the exact configuration of the
programming environment (such as internal solvers used by CVX).

For every parameter of the system, there is a different $\tau_{opt}$, with a
corresponding value of $\Gamma^{(L)}$ and $H^{(L)}_{LS}$ given by CVX. If we
want to explore a large parameter space of the system, it seems that we need a
different $\Gamma^{(L)}$ and $H^{(L)}_{LS}$ for each parameter point.
Surprisingly, we find that this is not always required. If
$\tau_{opt}\ll\delta$ for one set of parameters, we can substantially change
parameters of the system far from the qubits attached to baths, and still
obtain a value of $\tau\ll\delta$ with the same value of $\Gamma^{(L)}$ and
$H^{(L)}_{LS}$. This is shown in Fig. (6), where $\tau$ is calculated changing
various parameters away from the two qubits coupled with the bath, fixing
$H^{(L)}_{LS}=0$ and $\Gamma^{(L)}$ to be the same as in Table. 1. Over the
entire regime of chosen parameters $\tau\ll\delta$. Note that, in contrast to
previous plots, this is not be the optimal value of $\tau$. Nevertheless, if
$\tau\ll\delta$, we still get a completely positive Markovian QME preserving
local conservation laws and showing thermalization up to the chosen precision.
Given $\Gamma^{(L)}$ and $H^{(L)}_{LS}$, it is much easier to just check this
rather than finding the optimal value $\tau_{opt}$.

If parameters of the two qubits that are coupled to the bath are changed, we
can no longer use the same $\Gamma^{(L)}$ and $H^{(L)}_{LS}$. For example, if
we choose the same $\Gamma^{(L)}$ as in Fig. 6, and change $g_{1}$ to $0.2$
from $0.1$, we get $\tau=0.0014\gg\delta$.

The above observation suggests that the values of $\Gamma^{(L)}$ and
$H^{(L)}_{LS}$ obtained by CVX can be used to define a QME, independent of the
parameters in the bulk of the system. This is consistent with underlying
picture that each value of $\Gamma^{(L)}$ and $H^{(L)}_{LS}$ corresponds to a
different choice of the bath spectral function and the system-bath coupling
Hamiltonian. If we change any parameter of the qubits attached to the baths,
the change reflects substantially on the system-bath coupling Hamiltonian, so
the value of $\tau$ changes drastically from $\tau_{opt}$ obtained with
original parameters. If we change any parameter away from the two qubits, the
change reflects much less on the system-bath coupling Hamiltonian, causing
$\tau$ to be of the same order as the original value of $\tau_{opt}$. This
presents an exciting prospect for studying the dynamics of the system-bath
setup over a wide range of parameters, including a wide range of temperatures,
with physically consistent Markovian QMEs. Such studies may also be possible
for long chains, since local Markovian dissipation is favourable for tensor
network based numerical techniques. Such dissipation may, also, in principle,
be engineered in quantum computing and quantum simulation platforms, like ion
traps Barreiro  _et al._ (2011); Schindler  _et al._ (2013), Rydberg atoms
Weimer  _et al._ (2010); Nguyen  _et al._ (2018), superconducting qubits
García-Pérez  _et al._ (2020) and quantum dots Kim  _et al._ (2022).

Figure 6: $\tau$ vs $g_{4}$, for $N_{L}=2,N_{M}=4$, with
$\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$, $\beta=1$ and $g_{\ell}=0.1$ for
all $\ell$ unless otherwise mentioned. $\tau$ is computed from $\Gamma^{(L)}$
and $H^{(L)}_{LS}$ obtained from CVX for $N_{L}=2,N_{M}=4$, with
$\omega^{(\ell)}_{0}=1$, $\Delta_{\ell}=1$, $\beta=1$ and $g_{\ell}=0.1$ for
all $\ell$. The modified parameters for the plots are given by (i) (no
parameters changed), (ii) $\Delta_{3}=0.4,\Delta_{4}=1.2$, (iii)
$\omega^{(3)}_{0}=1.5,\omega^{(4)}_{0}=1.5,g_{5}=0.3$, (iv)
$\omega^{(3)}_{0}=1.5,\omega^{(4)}_{0}=1.5,g_{5}=0.3,\Delta_{4}=0.5$, (v)
$g_{3}=0.3$. We find that $\tau\ll\delta=10^{-6}$ even if parameters are
changed for qubits of the system that are not coupled to the baths.

##  V Summary and outlook

Searching for a physically consistent Markvian QME — A physically consistent
Markovian QME must satisfy complete positivity, obey local conservation laws
and be able to show thermalization. In this work, we have systematically gone
about searching for such QMEs. This is done in three steps, and the result in
each step has important consequences. Especially, we are led to introduce the
TOP problem, which is an optimization problem for finding a QME with all the
above properties up to a given precision. The TOP opens a completely new
avenue in the study of dissipative quantum systems.

We consider a finite-dimensional undriven system a part of which is weakly
coupled to a thermal bath. The microscopically derived QME written to leading
order in system-bath coupling is the RE, which is known to obey local
conservation laws and be able to show thermalization Tupkary  _et al._ (2022).
First, we show in generality that the RE violates complete positivity, unless
in extremely special cases. Although there are previous works showing this via
specific examples (for instance, Hartmann and Strunz (2020); Eastham  _et al._
(2016); Suárez  _et al._ (1992); Spohn (1980)), we are unaware of a model
independent proof similar to ours. Next, we prove that imposing complete
positivity and preservation of local conservation laws enforces the QME to be
of ‘local’ form. That is, the Lindblad operators and the Lamb-shift
Hamiltonian must have support only on the part of the system directly coupled
to the bath, and be identity elsewhere. This rules out the possibility of any
‘global’ forms of Lindblad equations, which are usually constructed to show
thermalization, to be consistent with local conservation laws. Then, we ask if
a ‘local’ Lindblad equation can be found which is able to show thermalization.
We find that, the task of finding such a Lindblad equation can be cast as an
optimization problem, which we call TOP. Most interestingly, this optimization
problem turns out to be a SDP. For given system and parameters, the SDP can be
efficiently solved using high-level programming packages like the CVX Matlab
package. The output of the TOP conclusively shows whether the desired type of
QME is possible for the chosen system parameters and temperature, up to a
chosen precision. For numerical example, we look at the TOP in a XXZ qubit
chain of few sites, fixing a reasonably high precision. When only the first
site is coupled to a bath, we find that, unless in extremes of temperatures,
it is impossible to find a local Lindblad equation that is capable of showing
thermalization up to the chosen precision.

Discussion in light of various existing forms of QMEs — Various forms of QMEs
have been derived it literature under various approximations (for example,
Becker  _et al._ (2022, 2021); Nathan and Rudner (2020); McCauley  _et al._
(2020); Harbola  _et al._ (2006); Nathan and Rudner (2022); Trushechkin
(2021); Trushechkin and Volovich (2016); Kleinherbers  _et al._ (2020);
Davidović (2020); Mozgunov and Lidar (2020); Kiršanskas  _et al._ (2018)),
along with the standard RE, local and eigenbasis Lindblad equations Breuer and
Petruccione (2006). Although the above example shows that there is no general
form of physically consistent Markovian QME, this does not immediately make
them unusable. Instead, it turns out that in each of these forms of QME, some
elements of the system density matrix are given correctly, while the others
are not Tupkary  _et al._ (2022). So, one needs to be careful in interpreting
the results from them, always keeping in mind their microscopic derivation and
approximations. The RE, despite not being completely positive, is provably
more accurate than all such Lindblad QMEs. To elucidate how this can happen,
imagine that, in a given setup, physically, the population of one energy
level, say, $\bra{E_{j}}\rho\ket{E_{j}}$, is zero in steady state. The RE
might then give a small negative value (say,
$\bra{E_{j}}\rho\ket{E_{j}}=-10^{-3}$), while any of the Lindblad equations
will give a larger positive value, which might be (say,
$\bra{E_{j}}\rho\ket{E_{j}}=0.1$). Either case is a problem if we want to
calculate various kinds of entropies, as often required in quantum information
and thermodynamics. In case of RE, unphysical results can often be ruled out
by checking the scaling with system-bath coupling Tupkary  _et al._ (2022);
Purkayastha  _et al._ (2020). This is often more difficult in Lindblad QMEs,
where approximations are often less controlled. The state obtained from the
recently derived ULE Nathan and Rudner (2020), which been shown to violate
local conservation laws Tupkary  _et al._ (2022), can be corrected to obtain
results as accurate as the RE Nathan and Rudner (2022). This re-instates the
local conservation laws, at the cost of also re-instating the same positivity
problem of the density matrix as in RE. In another recent work, a general form
of QME has been derived Becker  _et al._ (2022) which is more accurate than
RE, even though complete positivity of dynamics is still not guaranteed.

TOP and (non) Markovianity — In the microscopic picture, given the temperature
of the bath, the QME is completely defined by the bath spectral functions and
the type of system-bath coupling. The TOP can then be thought of as varying
over all possible bath spectral functions and types of system-bath couplings
to find the closest to satisfying thermalization the local Lindblad equation
can be. So, when TOP shows that the desired type of QME is impossible, it
means no matter what type of bath is attached and how it is coupled to the
system, for the chosen parameters, it is impossible to describe the dynamics
via a completely positive Markovian QME satisfying local conservation laws and
showing thermalization. The approach to thermal state must then have some non-
Markovian character for such system parameters and temperature. The output of
TOP, $\tau_{opt}$, shows non-trivial dependence on the system parameters and
the temperature. This dependence seems to capture how close to Markovian the
dynamics can be for the chosen parameters.

Surprises when two qubits are attached to bath — Surprisingly, we have found
that, when first two qubits of the few-site XXZ chain are attached to a bath,
solving the TOP shows that it is possible to find Lindbladians obeying local
conservation laws and showing thermalization up to quite high precision. This
holds over a considerable range of parameters, including a wide range of
temperatures. Notably, in this entire parameter regime, when one qubit was
coupled to a bath, such a QME was impossible.

Whenever the TOP shows a QME respecting all conditions is possible, standard
high-level programming packages used to solve the SDP also outputs one
possible form for such a QME. When two qubits are attached to the bath, the
form of QME so obtained, which respects all the requirements, is quite non-
trivial and would be hard to guess otherwise. Even more interestingly, we have
found that if we take one such QME obtained for one choice of system
parameters, and change some system parameters away from two qubits that couple
to the bath, the QME still satisfies all the requirements. This opens several
exciting possibilities that we describe below.

Future directions — Our results open the exciting possibility of studying the
dynamics of approach to thermal state in open quantum many-body systems using
physically consistent Markovian QMEs, over a wide range of parameters,
including a wide range of temperatures. This is particularly aided by the fact
that local Lindblad equations are favourable for tensor network techniques.
Studying such dynamics at finite temperatures is often quite challenging
otherwise, requiring simulation of non-Markovian dynamics Fux  _et al._
(2022); Purkayastha  _et al._ (2021); Lacerda  _et al._ (2022).

The TOP lets us find parameters of the system where local Lindblad equations
can show thermalization. For two qubits attached to bath, this range of
parameters can be considerably large, as we have seen. It may be possible to
design such local dissipation in quantum computing and quantum simulation
platforms like ion traps Barreiro  _et al._ (2011); Schindler  _et al._
(2013), Rydberg atoms Weimer  _et al._ (2010); Nguyen  _et al._ (2018),
superconducting qubits García-Pérez  _et al._ (2020) and quantum dots Kim  _et
al._ (2022). Especially in ion traps and Rydberg atom platforms, this offers
an interesting way to controllably prepare finite temperature states of
complex quantum many-body systems in these platforms, which is presently a
technological challenge. Usually, one would require global Lindblad
dissipators to ensure that a thermal state is prepared. This would be hard to
design in quantum simulation platforms if one wants to simulate complex many-
body systems. The possibility of having local dissipation confined to two
qubits offers a much easier alternative.

Moreover, as we have seen in the example of XXZ qubit chain, the dependence of
the output of the TOP, $\tau_{opt}$, on various parameters of the system
already encode rich and interesting physics. For complex quantum many-body
systems, one may need more scalable techniques for SDP, which is itself a
direction of research in computer science Majumdar  _et al._ (2020). Using
these techniques, the rich behavior of $\tau_{opt}$ with various parameters
can then be studied.

It is therefore clear that our results, especially the introduction of the
TOP, leads to new paradigm within the fields of quantum information,
computation and technology. Nevertheless, various questions still remain. One
main question concerns steady-state coherences Guarnieri  _et al._ (2018);
Purkayastha  _et al._ (2020); Trushechkin  _et al._ (2022); Cresser and Anders
(2021). When coupled to a thermal bath at any finite coupling, the system
density matrix will have coherences in energy eigenbasis of the system Cresser
and Anders (2021); Trushechkin  _et al._ (2022). These coherences can be
important in quantum information and thermodynamics Streltsov  _et al._
(2017); Francica  _et al._ (2019); Santos  _et al._ (2019); Narasimhachar and
Gour (2015) and are given correctly to the leading order by the RE Tupkary
_et al._ (2022); Purkayastha  _et al._ (2020); Cresser and Anders (2021).
However, it is not clear that the steady-state coherences calculated from
physically consistent Markovian QME obtained via TOP will be the same as those
obtained from RE. Further investigation is required in this respect, which
will be carried out in future works.

All code used in this work can be found at Tupkary (2022).

## 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). A.P
also acknowledges funding from the Danish National Research Foundation through
the Center of Excellence “CCQ” (Grant agreement no.: DNRF156).

##  Appendix A Casting Eq.(LABEL:RE6) to Eq.(15)

In this appendix we show the steps for taking Eq. (LABEL:RE6) to the form of
Eq. (15) which is more amenable to studying issues related to conservation of
complete positivity. We start with Eq. (LABEL:RE6), which we recall to be

$\displaystyle\frac{\partial{\rho}}{\partial
t}=i[\rho,H_{S}]-\epsilon^{2}\sum_{l}\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}\Big{\\{}a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}[F^{\dagger}_{\alpha},F_{\tilde{\alpha}}\rho]$
(46) $\displaystyle+c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}[\rho
F^{\dagger}_{\alpha},F_{\tilde{\alpha}}]+b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}[\rho
F_{\alpha}^{\dagger},F_{\tilde{\alpha}}]+a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}[F^{\dagger}_{\alpha},F_{\tilde{\alpha}}\rho]\Big{\\}}.$

This can be rewritten as

$\displaystyle\frac{\partial{\rho}}{\partial t}$
$\displaystyle=i[\rho,H_{S}]+\epsilon^{2}\sum_{l}\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}\Big{\\{}a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}\Big{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}-\frac{[F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho]}{2}\Big{)}+c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}\Big{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}+\frac{[F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho]}{2}\Big{)}$
(47)
$\displaystyle+b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}\Big{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}+\frac{[F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho]}{2}\Big{)}+a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}\Big{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}-\frac{[F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho]}{2}\Big{)}\Big{\\}},$

where, ${A,B}=AB+BA$ is the anti-commutator. Next, we note that the summation
$\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}$ in above equation, can be written
as
$\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}=\sum_{\alpha=\tilde{\alpha}=d^{2}}+\sum_{\alpha=1,\tilde{\alpha}=d^{2}}^{d^{2}-1}+\sum_{{\tilde{\alpha}}=1,\alpha=d^{2}}^{d^{2}-1}+\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}-1}.$
Using this, and the fact that $F_{d^{2}}=I_{S}/\sqrt{d}$ commutes with all
operators, we combine all commutator terms and write them as as
$i[\rho,H_{S}+H_{LS}]$ to obtain

$\displaystyle\frac{\partial{\rho}}{\partial t}$
$\displaystyle=i[\rho,H_{S}+H_{LS}]+\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}-1}\Gamma_{\alpha{\tilde{\alpha}}}\Big{(}F_{\tilde{\alpha}}\rho
F_{\alpha}^{\dagger}-\frac{\\{F_{\alpha}^{\dagger}F_{\tilde{\alpha}},\rho\\}}{2}\Big{)}.$
(48)

Here

$\displaystyle H_{LS}=\epsilon^{2}\sum_{l}\bigg{\\{}$ (49)
$\displaystyle\sum_{\alpha,{\tilde{\alpha}}=1}^{d^{2}}\Big{\\{}\frac{a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}}{2i}-\frac{c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}}{2i}-\frac{b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}}{2i}+\frac{a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}}{2i}\Big{\\}}F^{\dagger}_{\alpha}F_{\tilde{\alpha}}$
$\displaystyle+\sum_{\alpha=1}^{d^{2}-1}\frac{(a^{*}_{l\alpha}b_{l,d^{2}}+c^{\prime*}_{l\alpha}a^{\prime}_{ld^{2}}+b^{*}_{l\alpha}a_{l,d^{2}}+a^{\prime*}_{l\alpha}c^{\prime}_{l,d^{2}})}{2i\sqrt{d}}F^{\dagger}_{\alpha}$
$\displaystyle-\sum_{{\tilde{\alpha}}=1}^{d^{2}-1}\frac{(a^{*}_{ld^{2}}b_{l{\tilde{\alpha}}}+c^{\prime*}_{l,d^{2}}a^{\prime}_{l{\tilde{\alpha}}}+b^{*}_{l,d^{2}}a_{l{\tilde{\alpha}}}+a^{\prime*}_{l,d^{2}}c^{\prime}_{l{\tilde{\alpha}}})}{2i\sqrt{d}}F_{\tilde{\alpha}}\bigg{\\}}$

and

$\Gamma_{\alpha{\tilde{\alpha}}}=\epsilon^{2}\sum_{l}(a^{*}_{l\alpha}b_{l{\tilde{\alpha}}}+c^{\prime*}_{l\alpha}a^{\prime}_{l{\tilde{\alpha}}}+b^{*}_{l\alpha}a_{l{\tilde{\alpha}}}+a^{\prime*}_{l\alpha}c^{\prime}_{l{\tilde{\alpha}}}),$
(50)

$\alpha,\tilde{\alpha}$ going from $1$ to $d^{2}-1$. This is Eq. (15) given in
the main text.

##  Appendix B An example of RE violating complete positivity

In this section, we will present a simple example of the discussion in Sec.
II. Our setup consists of a two-qubit XXZ qubit chain, where only the first
qubit is connected to the bath modelled by an infinite number of bosonic
modes. Let $H$ be the Hamiltonian of the full set-up, given by

$H=H_{S}+\epsilon\,H_{SB}+H_{B},$ (51)

where

$\displaystyle H_{S}$
$\displaystyle=\frac{\omega_{0}}{2}(\sigma_{z}^{(1)}+\sigma_{z}^{(2)})$ (52)
$\displaystyle-g(\sigma_{x}^{(1)}\sigma_{x}^{(2)}+\sigma_{y}^{(1)}\sigma_{y}^{(2)}+\Delta\sigma_{z}^{(1)}\sigma_{z}^{(2)})$
$\displaystyle H_{SB}$
$\displaystyle=\sum_{r=1}^{\infty}(\kappa_{r}\hat{B}^{\dagger}_{r}\sigma^{(1)}_{-}+\kappa_{r}^{*}\hat{B}_{r}\sigma^{(1)}_{+})$
$\displaystyle H_{B}$
$\displaystyle=\sum_{r=1}^{\infty}\Omega_{r}\hat{B}^{\dagger}_{r}\hat{B}_{r}$

where $\sigma^{(\ell)}_{x,y,z}$ denotes the Pauli matrices acting on the
$\ell^{\text{th}}$ qubit,
$\sigma^{(\ell)}_{+}=(\sigma^{(\ell)}_{x}+i\sigma^{(\ell)}_{y})/2$,
$\sigma^{(\ell)}_{-}=(\sigma^{(\ell)}_{x}-i\sigma^{(\ell)}_{y})/2$,
$\hat{B}_{r}$ is bosonic annihilation operator for the $r^{\text{th}}$ mode of
the bath. Here, $\omega_{0}$, $g$, and $g\Delta$ represent the magnetic field,
the overall qubit-qubit coupling strength and the anisotropy respectively. The
RE for this setup can be computed to be Tupkary  _et al._ (2022)

$\displaystyle\frac{\partial{\rho}}{\partial
t}=i[\rho(t),H_{S}]+\epsilon^{2}\big{\\{}$
$\displaystyle[S^{\dagger},S^{(1)}\rho(t)]-[S^{\dagger},\rho(t)S^{(2)}]$ (53)
$\displaystyle+\text{H.c}\big{\\}}$

with

$\begin{aligned} S^{\dagger}&=\sigma_{+}^{(1)},\qquad S=\sigma_{-}^{(1)}\\\
S^{(1)}&=\sum_{j,k=1}^{4}\ket{E_{j}}\bra{E_{j}}\sigma^{(1)}_{-}\ket{E_{k}}\bra{E_{k}}D(j,k),\\\
S^{(2)}&=\sum_{j,k=1}^{4}\ket{E_{j}}\bra{E_{j}}\sigma^{(1)}_{-}\ket{E_{k}}\bra{E_{k}}C(j,k)\end{aligned},$
(54)

and

$\displaystyle C(j,k)$
$\displaystyle=\frac{\mathfrak{J}(E_{kj})n(E_{kj})}{2}-i\mathcal{P}\int_{0}^{\infty}d\omega\frac{\mathfrak{J}(\omega)n(\omega)}{\omega-
E_{kj}},$ (55) $\displaystyle D(j,k)$
$\displaystyle=\frac{e^{\beta(E_{kj}-\mu_{\ell})}\mathfrak{J}(E_{kj})n(E_{kj})}{2}$
$\displaystyle-i\mathcal{P}\int_{0}^{\infty}d\omega\frac{e^{\beta(\omega-\mu)}\mathfrak{J}(\omega)n(\omega)}{\omega-
E_{kj}},$ $\displaystyle\mathfrak{J}(\omega)$
$\displaystyle=\sum_{k=1}^{\infty}2\pi\left|\kappa_{k}\right|^{2}\delta(\omega-\Omega_{k}),$
$\displaystyle n(\omega)$ $\displaystyle=[e^{\beta\omega}-1]^{-1}.$

In above, $\mathfrak{J}(\omega)$ is called the bath spectral function. Let us
consider bosonic baths described by Ohmic spectral functions with Gaussian
cut-offs, $\mathfrak{J}(\omega)=\omega
e^{-(\omega/\omega_{c})^{2}}\Theta(\omega),$ where $\Theta(\omega)$ is the
Heaviside step function, and $\omega_{c}$ is the cut-off frequency. The above
operators can then be computed numerically.

The next step is to choose the basis $f_{i}$ and $g_{j}$ for operators on
$\mathcal{H}_{L}$ and $\mathcal{H}_{M}$. For the general case, one can start
with any set of linearly independent operators that forms a basis and includes
the identity operator, and then apply the Gram Schmidt orthonormalization
procedure to produce an orthonormal basis that includes the normalized
identity operator. For our case, one can easily verify that the set
$\\{-\sigma^{(i)}_{z}/\sqrt{2},\sigma^{(i)}_{-},\sigma^{(i)}_{+},I^{(i)}_{2}/\sqrt{2}\\}$
suffices, where $i=1$ for $\\{f_{i}\\}$ and $i=2$ for $\\{g_{j}\\}$, and
$I_{2}$ is the identity operator.

The basis for the full system $\\{F_{i}\\}$ can be constructed from the above
basis as described in subsection II.1.1, and is given by

$\displaystyle F_{i}$
$\displaystyle=f_{i}\otimes\frac{I_{2}^{(2)}}{2}\qquad\text{(for $i=1,2,3$)}$
(56) $\displaystyle F_{3i+j}$ $\displaystyle=f_{i}\otimes
g_{j}\qquad\text{(for $i=1,2,3,4$ and $j=1,2,3$)}$ $\displaystyle F_{16}$
$\displaystyle=\frac{I_{4}}{2}$

Any operator $X$ can be expanded in terms of the above basis as
$X=\sum_{\alpha}x_{\alpha}F_{\alpha}$, where
$x_{\alpha}=\braket{F_{\alpha},X}=\text{Tr}(F^{\dagger}_{\alpha}X)$. Thus,
expanding $S,S^{\dagger},S^{(1)},S^{(2)}$ [Eq. (54)], one can evaluate all the
coefficients in Eq. (13). Finally, one can compute the matrix $\Gamma$
according to Eq. (16). The matrix $\Gamma$ for this example, with parameters
chosen as $g=0.1,\omega_{0}=1,\omega_{c}=10,\beta=1,\mu=-0.5,\Delta=1$ is
given by

$\Gamma=\epsilon^{2}\setcounter{MaxMatrixCols}{15}\begin{bmatrix}0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&1.542+3.428i&0&0&0.014+0.047i&0&0&0&0&0&0&0&0&0\\\
0&1.542-3.428i&0&0&-0.18-0.007i&0&0.18+0.007i&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&-0.18+0.007i&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0.014-0.047i&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0.18-0.007i&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\end{bmatrix}$ (57)

We see that the above matrix has the expected structure of Eq. (17),

$\Gamma=\left[\begin{array}[]{ c | c }\Gamma_{\alpha,{\tilde{\alpha}}<4}&\Gamma_{\alpha<4,{\tilde{\alpha}}\geq 4}\\\ \hline\cr\Gamma_{\alpha\geq 4,{\tilde{\alpha}}<4}&0\\\ \end{array}\right]$ (58)

and crucially, $\Gamma_{\alpha<4,{\tilde{\alpha}}\geq 4}\neq 0$. Therefore, as
per the GKSL theorem, this RE will not preserve complete positivity. This
example was computed using QuTiP Johansson  _et al._ (2012, 2013).

##  Appendix C Effective Lindblad equation satisfying local conservation laws

In this appendix, we will show that Eq. (LABEL:cont_conjecture_2) implies Eq.
(31). The condition for a QME preserving complete positivity and obeying local
conservation laws is given by Eq. (LABEL:cont_conjecture_2), which we recall
to be

$\displaystyle-i[O_{M},H^{\prime}]+\sum_{\alpha_{M},{\tilde{\alpha}}_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}\Big{(}g_{\alpha_{M}}^{\dagger}O_{M}g_{{\tilde{\alpha}}_{M}}$
(59)
$\displaystyle-\frac{1}{2}O_{M}g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}-\frac{1}{2}g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}O_{M}\Big{)}=0,\qquad\forall
O_{M},$

where we write

$H^{\prime}=\sum_{\alpha_{M}=1}^{d_{M}^{2}}\nu_{d_{L}^{2},\alpha_{M}}g_{\alpha_{M}},$
(60)

for convenience. To move forward we will use the operator vector
correspondence from Ref. Watrous, 2018, where the vectorized version of the
operator $X$ is given by vec($X$), and can be constructed using linearity and

$\text{vec($\ket{i}\bra{j}$)}=\ket{i}\otimes\ket{j}^{*},$ (61)

where $\ket{j}^{*}$ denotes the complex conjugate of $\ket{j}$. We will apply
Eq. (61) to Eq. (LABEL:eq:cont_conjecture_3), using the identity (Eq. 1.132 of
Ref. Watrous, 2018)

$\text{vec($A_{0}BA_{1}^{T}$)}=(A_{0}\otimes A_{1})\text{vec($B$)}.$ (62)

Applying the vec operation on both sides of Eq. (LABEL:eq:cont_conjecture_3),
we obtain

$\displaystyle-i\,\text{vec($I_{M}O_{M}H^{\prime}$)}+i\,\text{vec($H^{\prime}O_{M}I_{M}$)}+\sum_{\alpha_{M},{\tilde{\alpha}}_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}$
(63)
$\displaystyle\Big{(}\text{vec($g_{\alpha_{M}}^{\dagger}O_{M}g_{{\tilde{\alpha}}_{M}}$)}-\frac{1}{2}\text{vec($I_{M}O_{M}g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}$)}$
$\displaystyle-\frac{1}{2}\text{vec($g_{\alpha_{M}}^{\dagger}g_{{\tilde{\alpha}}_{M}}O_{M}I_{M}$)}\Big{)}=0,\qquad\forall
O_{M}.$

Eq. (63) can be simplified using Eq. (62) to obtain,

$\displaystyle\bigg{\\{}-iI_{M}\otimes(H^{\prime})^{T}+iH^{\prime}\otimes
I_{M}\sum_{\alpha_{M},{\tilde{\alpha}}_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}$
(64) $\displaystyle\Big{(}g^{\dagger}_{\alpha_{M}}\otimes
g^{T}_{{\tilde{\alpha}}_{M}}-\frac{1}{2}I_{M}\otimes(g^{\dagger}_{\alpha_{M}}g_{{\tilde{\alpha}}_{M}})^{T}$
$\displaystyle-\frac{1}{2}(g^{\dagger}_{\alpha_{M}}g_{{\tilde{\alpha}}_{M}})\otimes
I_{M}\Big{)}\bigg{\\}}\text{vec($O_{M}$)}=0,\qquad\forall O_{M}.$

Eq. (64) is of the form $\mathcal{M}\,\text{vec($O_{M}$)}=0$ for all hermitian
$O_{M}$. Since hermitian matrices (such as $O_{M}$) form a basis for the
entire space of operators, this implies $\mathcal{M}\,\text{vec($X$)}=0$ for
all operators $X$. This is because one can expand $X$ as a linear combination
of hermitian operators (such as $O_{M}$). Now,
$\mathcal{M}\,\text{vec($X$)}=0\,\forall\,X$ implies $\mathcal{M}=0$.
Therefore, Eq. (64) implies

$\displaystyle\mathcal{M}=$ $\displaystyle-
iI_{M}\otimes(H^{\prime})^{T}+iH^{\prime}\otimes
I_{M}+\sum_{\alpha_{M},{\tilde{\alpha}}_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},{\tilde{\alpha}}_{M}}$
(65) $\displaystyle\Big{(}g^{\dagger}_{\alpha_{M}}\otimes
g^{T}_{{\tilde{\alpha}}_{M}}-\frac{1}{2}I_{M}\otimes(g^{\dagger}_{\alpha_{M}}g_{{\tilde{\alpha}}_{M}})^{T}$
$\displaystyle-\frac{1}{2}(g^{\dagger}_{\alpha_{M}}g_{{\tilde{\alpha}}_{M}})\otimes
I_{M}\Big{)}=0.$

If $\mathcal{M}=0$, then $\text{Tr}(\mathcal{M})=0$. Taking the trace of Eq.
(65), and using the orthonormality of $\\{g_{i}\\}$ along with the fact that
$\text{Tr}(g_{i})=\delta_{i,d^{2}_{M}}$, we obtain

$\text{Tr}(\mathcal{M})=-\sum_{\alpha_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},\alpha_{M}}d_{M}=0$
(66)

which implies

$\displaystyle\sum_{\alpha_{M}=1}^{d_{M}^{2}-1}\widetilde{\Lambda}_{\alpha_{M},\alpha_{M}}=\sum_{\alpha_{M}=1}^{d_{M}^{2}-1}\sum_{\alpha_{L}=1}^{d_{L}^{2}}\widetilde{\Gamma}_{(\alpha_{L},\alpha_{M}),(\alpha_{L},\alpha_{M})}=0.$
(67)

which is Eq. (31) in the main text.

##  Appendix D The condition for thermalization

From fundamental principles of quantum statistical mechanics, we expect the
system to thermalize when coupled to baths at equal temperatures. The exact
condition that QME’s must obey to satisfy thermalization has been derived in
Ref. Tupkary  _et al._ , 2022. For the sake of completeness, we recall that
discussion here.

Let the total system Hamiltonian be given by $H=H_{S}+\epsilon H_{SB}+H_{B}$,
where $\epsilon$ is a dimensionless parameter controlling the strength of the
system bath coupling, and $H_{SB}$ is the system bath coupling Hamiltonian. We
proceed by obtaining an order-by-order solution to the steady state of our
QME. Any QME describing our setup can be expanded in the so-called time-
convolution-less form Breuer and Petruccione (2006),

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

where $\mathcal{L}$ could in general be time-dependent operators and
$L_{0}(t)[\rho(t)]=i[\rho(t),H_{S}]$. For quantum master equations written to
second-order in system-bath coupling, the above summation can be truncated at
second order. Denoting
$\mathcal{L}_{2m}\equiv\lim_{t\rightarrow\infty}\mathcal{L}_{2m}(t)$, the
steady state $\rho_{SS}$ can be given by

$\rho_{SS}=\lim_{t\to\infty}e^{t(\mathcal{L}_{0}+\epsilon^{2}\mathcal{L}_{2})}\rho(0),$
(69)

which is assumed to be unique. The steady state satisfies

$0=\sum_{m=0}^{\infty}\epsilon^{2m}\mathcal{L}_{2m}[\rho_{SS}].$ (70)

We can then perform an expansion of $\rho_{SS}$ in the even powers of
$\epsilon$ as

$\rho_{SS}=\sum_{m=0}^{\infty}\epsilon^{2m}\rho^{(2m)}_{SS}$ (71)

Using Eq. (71) in Eq. (70), we can obtain an order by order solution of
$\rho_{SS}$. At the zeroth order in $\epsilon$, we obtain

$[\rho^{(0)}_{SS},H_{S}]=0.$ (72)

Assuming that the Hamiltonian has no degeneracies, Eq. (72) implies that
$\rho^{(0)}_{SS}$ is diagonal in the energy eigenbasis,

$\rho^{(0)}_{SS}=\sum_{i}p_{i}\ket{E_{i}}\bra{E_{i}}.$ (73)

where $\ket{E_{i}}$ is an eigenstate of the system. At second order in
$\epsilon$ ($m=1$), we obtain the following two equations,

$\bra{E_{i}}\mathcal{L}_{2}[\rho^{(0)}_{SS}]\ket{E_{i}}=0,\quad\forall i$ (74)
$\displaystyle i(E_{i}-E_{j})\bra{E_{i}}\rho^{(2)}_{SS}\ket{E_{j}}$ (75)
$\displaystyle+\epsilon^{2}\bra{E_{i}}\mathcal{L}_{2}[\rho^{(0)}_{SS}]\ket{E_{j}}$
$\displaystyle=0,\quad\forall i\neq j$

Since $\rho^{(0)}_{SS}$ is diagonal in the energy eigenbasis, Eq. (74)
determines the diagonal elements of $\rho^{(0)}_{SS}$. Having obtained
$\rho^{(0)}_{SS}$, Eq. (75) then determines the off-diagonal elements of
$\rho^{(2)}_{SS}$. Note from above equations that the leading order diagonal
elements of $\rho_{SS}$ are independent of $\epsilon$. It can also be shown
that the leading order off-diagonal elements of $\rho_{SS}$ in the energy
eigenbasis of the system scale as $\epsilon^{2}$. As discussed in the main
text, the QME thermalizes if

$\lim_{\epsilon\to 0}\rho_{SS}=\rho_{th}$ (76)

where $\rho_{th}$ is the Gibbs state of the system given by

$\rho_{th}=\frac{e^{-\beta H_{S}}}{\text{Tr}[e^{-\beta H_{S}}]}.$ (77)

We then conclude that the thermalization in this sense is a statement about
leading order diagonal elements of $\rho_{SS}$. Substituting Eq. (76) in Eq.
(74), we obtain the following condition on $\mathcal{L}_{2}$ for the system to
thermalize,

$\bra{E_{i}}\mathcal{L}_{2}[\rho_{th}]\ket{E_{i}}=0\quad\forall i.$ (78)

##  Appendix E Semidefinite Programming (SDP)

###  E.1  Basic Theory

In this section, we present the theoretical framework of semidefinite
programming (SDP). We follow the definition of SDPs given in page 57 of Ref.
Watrous, 2018. In what follows, we will use $\Phi$ and $\Psi$ to denote
hermitian preserving linear maps. We will also use $\Phi^{\dagger}$ to denote
the “adjoint map” Watrous (2018), which is defined as the unique linear map
that satisfies

$\braket{A,\Phi(B)}=\braket{\Phi^{\dagger}(A),B}$ (79)

where

$\braket{A,B}=\text{Tr}(A^{\dagger}B)$ (80)

denotes the Hilbert Schmidt inner product.

An SDP is defined by the tuple ($\Phi$, $\Psi$, $A$, $B$, $C$), where
$\Phi,\Psi$ are hermitian-preserving linear maps, and $A,B,C$ are hermitian
operators. The “primal” problem of the SDP is given by

$\displaystyle\text{maximize}:$ $\displaystyle\braket{A,X}\quad w.r.t.\hskip
5.0ptX$ (81) $\displaystyle\text{subject to}:$
$\displaystyle\Phi(X)=B,~{}\Psi(X)\leq C,~{}X\geq 0,$

where the inequalities represent matrix inequalities. I.e, $A\geq B$ is
equivalent to $A-B\geq 0$ and implies $A-B$ is positive semidefinite. We will
use the notation $X_{f}$ to denote any “feasible” value of $X$ that satisfies
the three constraints in Eq. (81), and $\mathcal{P}$ to denote the maximum
value of $\braket{A,X}$ attained in Eq. (81) (assuming there is atleast one
$X$ which satisfies constraints).

For every “primal” problem, there exists a “dual” problem given by

$\displaystyle\text{minimize}:$ $\displaystyle\braket{B,Y}+\braket{C,Z}\quad
w.r.t.\hskip 5.0ptY,Z$ (82) $\displaystyle\text{subject to}:$
$\displaystyle\Phi^{\dagger}(Y)+\Psi^{\dagger}(Z)\geq A,$ $\displaystyle
Y\text{is hermitian},Z\geq 0.$

We will use the notation $(Y_{f},Z_{f})$ to denote any “feasible” value of
$Y,Z$ that satisfies the three constraints in Eq. (82), and $\mathcal{D}$ to
denote the minimum value of $\braket{B,Y}+\braket{C,Z}$ attained in Eqs. (82)
(assuming atleast some $(Y,Z)$ satisfies constraints).

Semidefinite programs have a notion of duality associated with them, which
relates properties of the primal and the dual problems. In particular, it can
be shown that

$\mathcal{P}\leq\mathcal{D}$ (83)

a property known as “weak duality”. In most situations, it can be shown that
$\mathcal{P}=\mathcal{D}$, i.e, equality holds in Eq. (83). This condition is
known as “strong duality”.

By weak duality and the definition of our primal and dual problems, using Eq.
(81),(82), and (83), we obtain

$\braket{A,X_{f}}\leq\mathcal{P}\leq\mathcal{D}\leq\braket{B,Y_{f}}+\braket{C,Z_{f}}.$
(84)

From Eq. (84), any feasible choice of inputs to the primal and dual problem
($X_{f},Y_{f},Z_{f}$) leads to lower and upper bounds on the optimal values of
the primal and dual problems [see Fig. (7)]. In particular, if we obtain
$\braket{A,X_{f}}=\braket{B,Y_{f}}+\braket{C,Z_{f}}$, equality holds
throughout in Eq. (84). This property can therefore be exploited to obtain
exact solutions to the primal problem of an SDP.

Figure 7: Schematic representing weak duality for the SDP given in Eqs. (81)
and (82), according to Eq. (84). $X_{f}^{(j)}$ and $(Y_{f}^{(j)},Z_{f}^{(j)})$
represents any feasible input to the primal and dual problems respectively.
Any such inputs yield upper and lower bounds on the solutions of the primal
and dual problems.

We will show in Sec. E.2 that the thermal optimization problem (TOP) in Eq.
(40) can be reduced to the standard form of SDP [Eq. (81)].

###  E.2  Reducing the thermal optimization problem (TOP) to standard form

Recall that the TOP was given by [Eq.(40)]

minimize : $\displaystyle\tau$ (85) subject to : $\displaystyle
H^{(L)}_{LS}\text{ is
hermitian,}~{}\text{Tr}(\Gamma^{(L)})=1,~{}\Gamma^{(L)}\geq 0.$

See Eq. (39) for definition of $\tau$ and Eq.
(LABEL:eq:local_lindblad_for_TOP) for definition of $H^{(L)}_{LS}$,
$\Gamma^{(L)}$. In this subsection, we will show how the TOP from Eq. (85) can
be reduced to the standard form of an SDP. We note that the standard form of
SDP in Eq. (81) is not yet suitable for this purpose. Therefore, we replace
$A\rightarrow-A,C\rightarrow-C,\Psi\rightarrow-\Psi,Y\rightarrow-Y,$ in Eq.
(81) and Eq. (82), leaving $\Phi$, $B$, $X$ and $Z$ unchanged. Since
maximizing any function is the same as minimizing its negative, we obtain the
new “primal” form as

$\displaystyle\text{minimize}:$ $\displaystyle\braket{A,X}$ (86)
$\displaystyle\text{subject to}:$ $\displaystyle\Phi(X)=B,~{}\Psi(X)\geq
C,~{}X\geq 0,$

where we use $\widetilde{P}$ to denote the minimum value of $\braket{A,X}$
obtained in Eq. (86). The new “dual” form is written as

$\displaystyle\text{maximize}:$ $\displaystyle\braket{B,Y}+\braket{C,Z}$ (87)
$\displaystyle\text{subject to}:$
$\displaystyle\Phi^{\dagger}(Y)+\Psi^{\dagger}(Z)\leq A,$ $\displaystyle
Y\text{is hermitian},~{}Z\geq 0,$

where we use $\widetilde{D}$ to denote the minimum value of
$\braket{B,Y}+\braket{C,Z}$ obtained in Eq. (87). Eq. (84) is then transformed
into [see Fig (8)],

$\braket{A,X_{f}}\geq\widetilde{\mathcal{P}}\geq\widetilde{\mathcal{D}}\geq\braket{B,Y_{f}}+\braket{C,Z_{f}}.$
(88)

Figure 8: Schematic representing weak duality for the SDP given in Eqs. (86)
and (87), according to Eq. (88). $X_{f}^{(j)}$ and $(Y_{f}^{(j)},Z_{f}^{(j)})$
represents any feasible input to the primal and dual problems respectively.
Any such inputs yield upper and lower bounds on the solutions of the primal
and dual problems.

We will now show how to reduce the TOP from Eq. (85) to Eq. (86), via a series
of changes to the optimization problem in Eq. (86). We do so in three steps.

Step 1: Our first step is to write down a primal optimization problem whose
solution (minimum value attained) is equal to

$||K||_{1}=\text{Tr}(\sqrt{K^{\dagger}K}).$ (89)

for any hermitian matrix $K$. Let $\Pi_{p}$ and $\Pi_{n}$ be projectors onto
the positive and negative eigenspaces of $K$. In this case,

$||K||_{1}=\text{Tr}(\Pi_{p}K\Pi_{p})-\text{Tr}(\Pi_{n}K\Pi_{n}).$ (90)

Let us now consider the optimization problem given by,

$\displaystyle\text{minimize}:$
$\displaystyle\left\langle\begin{pmatrix}I&0\\\
0&I\end{pmatrix},\begin{pmatrix}P&.\\\ .&Q\end{pmatrix}\right\rangle$ (91)
$\displaystyle\text{subject to}:$ $\displaystyle\Psi_{1}\begin{pmatrix}P&.\\\
.&Q\end{pmatrix}=\begin{pmatrix}P&0\\\
0&Q\end{pmatrix}\geq\begin{pmatrix}K&0\\\ 0&-K\end{pmatrix},$
$\displaystyle\begin{pmatrix}P&.\\\ .&Q\end{pmatrix}\geq 0,$

where we use dots to represent arbitrary blocks of the matrices which can
always be set to zero without affecting the objective function or constraints.
Note that $\Psi_{1}$ in Eq. (91) is a map that replaces the off-diagonal
blocks with null matrices, leaving the diagonal blocks unchanged. It is easy
to see that Eq. (91) is of the form Eq. (86) (with $\Phi$ and $B$ omitted
i.e., no equality constraint). Thus Eq. (91) is an SDP.

The dual problem to the primal problem in Eq. (91) is given by

$\displaystyle\text{maximize}:$
$\displaystyle\left\langle\begin{pmatrix}K&0\\\
0&-K\end{pmatrix},\begin{pmatrix}\bar{P}&.\\\
.&\bar{Q}\end{pmatrix}\right\rangle$ (92) $\displaystyle\text{subject to}:$
$\displaystyle\Psi_{1}^{\dagger}\begin{pmatrix}\bar{P}&.\\\
.&\bar{Q}\end{pmatrix}=\begin{pmatrix}\bar{P}&0\\\
0&\bar{Q}\end{pmatrix}\leq\begin{pmatrix}I&0\\\ 0&I\end{pmatrix}$
$\displaystyle\begin{pmatrix}\bar{P}&.\\\ .&\bar{Q}\end{pmatrix}\geq 0.$

where $\Psi_{1}^{\dagger}$ turns out to be the same map as $\Psi_{1}$ using
Eq. (79). It can be seen that Eq. (92) is of the form Eq. (87) (again with
$\Phi$ and $B$ omitted).

We will now show that the optimal values attained in both the primal and dual
problems in Eqs. (91) and (92) is equal to $||K||_{1}$. To show this, note
that setting

$P_{f}=\Pi_{p}K\Pi_{p},\quad Q_{f}=-\Pi_{n}K\Pi_{n}$ (93)

(where $P_{f}$ and $Q_{f}$ denote ‘feasible’ choices of $P$ and $Q$
respectively) in Eq. (91) allows us to obtain $||K||_{1}$ in the primal
objective function. Furthermore, setting

$\bar{P}_{f}=\Pi_{p},\quad\bar{Q}_{f}=\Pi_{n}$ (94)

in Eq. (92) allows us to obtain $||K||_{1}$ in the dual objective function.
Thus, we have explicitly constructed feasible choices of inputs to the primal
and dual problems of Eqs. (91) and (92) respectively, that yield $||K||_{1}$
in the primal and dual objective functions. Therefore, according to Eq. (88),
the optimal values attained in the primal and dual problems are both exactly
equal to $||K||_{1}$.

Step 2: In the first step we constructed an SDP whose solution is equal to
$||K||_{1}$, for a fixed $K$. We will now construct an SDP which computes the
minimum value of $||K||_{1}$, subject to some constraints on $K$. We will
first recast the problem in Eq. (91) as

$\displaystyle\text{minimize}:$ $\displaystyle\text{Tr}(P)+\text{Tr}(Q)$ (95)
$\displaystyle\text{subject to}:$ $\displaystyle P\geq K,~{}Q\geq-K,~{}P,Q\geq
0.$

Eq. (91) computes

$||K||_{1}=\sum_{i}|K_{ii}|.$ (96)

when $K$ is diagonal. Let $\mathcal{G}$ be a linear, hermitian preserving map
that always outputs a diagonal matrix. Then, the optimization problem given by

$\displaystyle\text{minimize}:$ $\displaystyle\text{Tr}(P)+\text{Tr}(Q)$ (97)
$\displaystyle\text{subject to}:$ $\displaystyle
P\geq\mathcal{G}(R),~{}Q\geq-\mathcal{G}(R),$
$\displaystyle\Phi(R)=B,~{}~{}P,Q,R\geq 0,$

computes $\min_{R\geq 0,\Phi(R)=B}\sum_{i}|\mathcal{G}(R)_{ii}|$. We will now
begin to connect the above formalism to TOP from Eq. (85). We will show how
$R$ can be chosen to reflect the optimization over $H^{(L)}_{LS}$ and
$\Gamma^{(L)}$, and $\Phi$ can be chosen to reflect the trace constraint on
$\Gamma^{(L)}$. Then, we will specify a map $\mathcal{G}$ that takes $R$ as
input (i.e, $H^{(L)}_{LS}$ and $\Gamma^{(L)}_{LS}$), and outputs a diagonal
matrix such that the objective function computes

$\tau=\sum_{i}\left|\bra{E_{i}}\mathcal{L}_{2}(\rho_{th})\ket{E_{i}}\right|.$
(98)

Step 3: In the thermal optimization problem [Eq. (85)], we have an
optimization over $\Gamma^{(L)}\geq 0$, and hermitian $H^{(L)}_{LS}$. We use
the fact that any hermitian matrix $H^{(L)}_{LS}$ can be written as a
$H^{(L)}_{LS}=S-T$, where $S,T\geq 0$. Moreover $S-T$ for any $S,T\geq 0$ is
always hermitian. We will now replace the hermitian $H^{(L)}_{LS}$ with the
difference of positive matrices $S-T$. This is needed, since semidefinite
programs can only handle optimization over positive semidefinite variables.
Furthermore, let us identify $\Gamma^{(L)}$ with some matrix $U$. Now consider
the map $\mathcal{G}$ that acts as follows :

$\mathcal{G}\begin{pmatrix}S&.&.\\\ .&T&.\\\
.&.&U\end{pmatrix}\equiv\mathcal{G}(S,T,U)\equiv\begin{pmatrix}\bra{E_{1}}\mathcal{L}_{2}(\rho_{th})\ket{E_{1}}&0&\ldots&0\\\
0&\bra{E_{2}}\mathcal{L}_{2}(\rho_{th})\ket{E_{2}}&\ldots&0\\\
\vdots&\ldots&\ddots&\vdots\\\
0&\ldots&0&\bra{E_{d}}\mathcal{L}_{2}(\rho_{th})\ket{E_{d}},\end{pmatrix}$
(99)

where the map constructs $\mathcal{L}_{2}$ according to Eq.
(LABEL:eq:local_lindblad_for_TOP) after setting $H^{(L)}_{LS}=S-T$, and
$\Gamma^{(L)}=U$. Now, we consider a specific case of the optimization problem
in Eq. (97), for the choice of $\mathcal{G}$ in Eq. (99). We obtain,

$\displaystyle\text{minimize}:$ $\displaystyle\text{Tr}(P)+\text{Tr}(Q)$ (100)
$\displaystyle\text{subject to}:$ $\displaystyle
P\geq\mathcal{G}\begin{pmatrix}S&.&.\\\ .&T&.\\\
.&.&U\end{pmatrix},~{}Q\geq-\mathcal{G}\begin{pmatrix}S&.&.\\\ .&T&.\\\
.&.&U\end{pmatrix},$ $\displaystyle\Phi_{1}\begin{pmatrix}S&.&.\\\ .&T&.\\\
.&.&U\end{pmatrix}$ $\displaystyle=\text{Tr}(U)=1,P,Q,S,T,U\geq 0.$

Since $\mathcal{G}$ always outputs a diagonal matrix [see Eq. (99)], Eq. (100)
computes

$\displaystyle\min$
$\displaystyle\sum_{i}\left|\mathcal{G}\begin{pmatrix}S&&\\\ &T&\\\
&&U\end{pmatrix}_{ii}\right|$ (101) subject to : $\displaystyle S,T,U\geq
0,\text{Tr}(U)=1$

Since $H^{(L)}_{LS}$ can always be written as $S-T$, and $\Gamma^{(L)}$ as
$U$, Eq. (101) [and therefore Eq. (100) ] is identical to the thermal
optimization problem in Eq. (85). Therefore, Eq. (100) computes $\tau_{opt}$.

All that remains is converting Eq. (100) to the standard form Eq. (86). Eq.
(100) can be obtained from Eq. (100) after choosing

$\displaystyle A=\begin{pmatrix}I&&&&\\\ &I&&&\\\ &&0&&\\\ &&&0&\\\ &&&&0\\\
\end{pmatrix},\quad B=1,\quad X=\begin{pmatrix}P&.&.&.&.\\\ .&Q&.&.&.\\\
.&.&S&.&.\\\ .&.&.&T&.\\\ .&.&.&.&U\\\ \end{pmatrix},\quad C=0,$ (102)
$\displaystyle\Psi\begin{pmatrix}P&.&.&.&.\\\ .&Q&.&.&.\\\ .&.&S&.&.\\\
.&.&.&T&.\\\ .&.&.&.&U\\\
\end{pmatrix}=\begin{pmatrix}P-\mathcal{G}(S,T,U)&0\\\
0&Q+\mathcal{G}(S,T,U)\\\ \end{pmatrix},\quad\Phi\begin{pmatrix}P&.&.&.&.\\\
.&Q&.&.&.\\\ .&.&S&.&.\\\ .&.&.&T&.\\\ .&.&.&.&U\\\
\end{pmatrix}=\text{Tr}(U)$

Recall that it is helpful to think of $P,Q$ as variables needed to compute the
objective function $\tau$ [Eq. (85)], $S,T$ are variables that give rise to
$H^{(L)}_{LS}=S-T$, and $U$ is a variable that encodes $\Gamma^{(L)}$. We have
therefore shown that the TOP [Eq. (85)] is an SDP.

It is to be noted that it is not necessary to reduce the TOP to the standard
form of an SDP in order to use CVX Grant and Boyd (2014). Infact, the TOP from
Eq. (85) can be directly implemented into CVX, which itself handles the
construction of the dual problem automatically in the background.

## References

  * 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).
  * Hur  _et al._ (2016) K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró, Comptes Rendus Physique 17, 808 (2016).
  * 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).
  * Breuer and Petruccione (2006) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2006).
  * Rivas and Huelga (2012) Á. Rivas and S. F. Huelga, SpringerBriefs in Physics (2012).
  * Carmichael (2002) H. Carmichael,  _Statistical Methods in Quantum Optics 1. Master Equations and Fokker–Planck Equations_ (Springer-Verlag Berlin Heidelberg, 2002).
  * Watrous (2018) J. Watrous,  _The Theory of Quantum Information_ (Cambridge University Press, 2018).
  * Siddhu and Tayur (2022) V. Siddhu and S. Tayur, “Five starter pieces: Quantum information science via semidefinite programs,” in  _Tutorials in Operations Research: Emerging and Impactful Topics in Operations_ (2022) Chap. 3, pp. 59–92.
  * Mazziotti (2004) D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
  * Mazziotti (2011) D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
  * Wolf  _et al._ (2008) M. M. Wolf, J. Eisert, T. S. Cubitt, and J. I. Cirac, Phys. Rev. Lett. 101, 150402 (2008).
  * Cubitt  _et al._ (2012) T. S. Cubitt, J. Eisert, and M. M. Wolf, Communications in Mathematical Physics 310, 383 (2012).
  * Kiilerich and Mølmer (2018) A. H. Kiilerich and K. Mølmer, Phys. Rev. A 97, 052113 (2018).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Gorini  _et al._ (1978) V. Gorini, A. Frigerio, M. Verri, A. Kossakowski, and E. Sudarshan, Reports on Mathematical Physics 13, 149 (1978).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * 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).
  * Tupkary  _et al._ (2022) D. Tupkary, A. Dhar, M. Kulkarni, and A. Purkayastha, Phys. Rev. A 105, 032208 (2022).
  * 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).
  * Nathan and Rudner (2020) F. Nathan and M. S. Rudner, Phys. Rev. B 102, 115109 (2020).
  * Spohn (1980) H. Spohn, Rev. Mod. Phys. 52, 569 (1980).
  * Grant and Boyd (2014) M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx (2014).
  * Walls and Milburn (2008) D. Walls and G. J. Milburn,  _Quantum optics_ (Springer-Verlag Berlin Heidelberg, 2008).
  * Becker  _et al._ (2021) T. Becker, L.-N. Wu, and A. Eckardt, Phys. Rev. E 104, 014110 (2021).
  * Prosen and Žnidarič (2009) T. Prosen and M. Žnidarič, Journal of Statistical Mechanics: Theory and Experiment 2009, P02035 (2009).
  * Mendoza-Arenas  _et al._ (2015) J. J. Mendoza-Arenas, S. R. Clark, and D. Jaksch, Phys. Rev. E 91, 042129 (2015).
  * Žnidarič  _et al._ (2010) M. Žnidarič, T. c. v. Prosen, G. Benenti, G. Casati, and D. Rossini, Phys. Rev. E 81, 051135 (2010).
  * Barreiro  _et al._ (2011) J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 (2011).
  * Schindler  _et al._ (2013) P. Schindler, M. Müller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, Nature Physics 9, 361 (2013).
  * Weimer  _et al._ (2010) H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Nguyen  _et al._ (2018) T. L. Nguyen, J. M. Raimond, C. Sayrin, R. Cortiñas, T. Cantat-Moltrecht, F. Assemat, I. Dotsenko, S. Gleyzes, S. Haroche, G. Roux, T. Jolicoeur, and M. Brune, Phys. Rev. X 8, 011032 (2018).
  * García-Pérez  _et al._ (2020) G. García-Pérez, M. A. C. Rossi, and S. Maniscalco, npj Quantum Information 6, 1 (2020).
  * Kim  _et al._ (2022) C. W. Kim, J. M. Nichol, A. N. Jordan, and I. Franco, PRX Quantum 3, 040308 (2022).
  * Becker  _et al._ (2022) T. Becker, A. Schnell, and J. Thingna, Phys. Rev. Lett. 129, 200403 (2022).
  * McCauley  _et al._ (2020) G. McCauley, B. Cruikshank, D. I. Bondar, and K. Jacobs, npj Quantum Information 6 (2020).
  * Harbola  _et al._ (2006) U. Harbola, M. Esposito, and S. Mukamel, Phys. Rev. B 74, 235309 (2006).
  * Nathan and Rudner (2022) F. Nathan and M. S. Rudner,  (2022), 10.48550/ARXIV.2206.02917.
  * Trushechkin (2021) A. Trushechkin, Physical Review A 103 (2021), 10.1103/physreva.103.062226.
  * Trushechkin and Volovich (2016) A. S. Trushechkin and I. V. Volovich, EPL (Europhysics Letters) 113, 30005 (2016).
  * 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).
  * Kiršanskas  _et al._ (2018) G. Kiršanskas, M. Franckié, and A. Wacker, Phys. Rev. B 97, 035432 (2018).
  * Purkayastha  _et al._ (2020) A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, and J. Goold, npj Quantum Information 6 (2020).
  * Fux  _et al._ (2022) G. E. Fux, D. Kilda, B. W. Lovett, and J. Keeling,  (2022), 10.48550/ARXIV.2201.05529.
  * Purkayastha  _et al._ (2021) A. Purkayastha, G. Guarnieri, S. Campbell, J. Prior, and J. Goold, Phys. Rev. B 104, 045417 (2021).
  * Lacerda  _et al._ (2022) A. M. Lacerda, A. Purkayastha, M. Kewming, G. T. Landi, and J. Goold,  (2022), 10.48550/ARXIV.2206.01090.
  * Majumdar  _et al._ (2020) A. Majumdar, G. Hall, and A. A. Ahmadi, Annual Review of Control, Robotics, and Autonomous Systems 3, 331 (2020).
  * Guarnieri  _et al._ (2018) G. Guarnieri, M. Kolář, and R. Filip, Phys. Rev. Lett. 121, 070401 (2018).
  * Trushechkin  _et al._ (2022) A. S. Trushechkin, M. Merkli, J. D. Cresser, and J. Anders, AVS Quantum Science 4, 012301 (2022).
  * Cresser and Anders (2021) J. D. Cresser and J. Anders, Phys. Rev. Lett. 127, 250601 (2021).
  * Streltsov  _et al._ (2017) A. Streltsov, G. Adesso, and M. B. Plenio, Rev. Mod. Phys. 89, 041003 (2017).
  * Francica  _et al._ (2019) G. Francica, J. Goold, and F. Plastina, Phys. Rev. E 99, 042105 (2019).
  * Santos  _et al._ (2019) J. P. Santos, L. C. Céleri, G. T. Landi, and M. Paternostro, npj Quantum Information 5 (2019).
  * Narasimhachar and Gour (2015) V. Narasimhachar and G. Gour, Nature Communications 6 (2015).
  * Tupkary (2022) D. Tupkary, https://github.com/dtupkary/SearchingLindbladians (2022).
  * 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).