# An optimal superconducting hybrid machine

Rosa López<EMAIL_ADDRESS>Instituto de Física Interdisciplinar y
Sistemas Complejos IFISC (CSIC-UIB), E-07122 Palma de Mallorca, Spain  Jong
Soo Lim<EMAIL_ADDRESS>Department of molecular design, Arontier Co.,
15F Sewon Bldg., Gangnam-daero 241, Seocho-gu, Seoul, Republic of Korea  Kun
Woo Kim<EMAIL_ADDRESS>Department of Physics, Chung-Ang University, 06974
Seoul, Republic of Korea

###### Abstract

Optimal engine performances are accomplished by quantum effects. Here we
explore two routes towards ideal engines, namely (1) quantum systems that
operate as hybrid machines being able to perform more than one useful task and
(2) the suppression of fluctuations in doing such tasks. For classical
devices, the absence of fluctuations is conditioned by a high entropy
production as dictate the thermodynamic uncertainty relations. Here we
generalize such relations for multiterminal conductors that operate as hybrid
thermal machines. These relations are overcome in quantum conductors as we
demonstrate for a double quantum dot contacted to normal metals and a
reservoir being a generator of entangled Cooper pairs.

_Introduction—_ The second law of thermodynamics dictates that Carnot
efficiency is the maximum efficiency for a thermal machine that delivers zero
power and works reversibly [1]. From a practical point of view, however, under
a nonequilibrium situation a thermal engine generates some power and a finite
amount of entropy. In such a scenario, a new kind of nonequilibrium principle
establishes a connection between power fluctuations and entropy production,
which is called the thermodynamic uncertainty relation (TUR) [2, 3, 4, 5, 6,
7, 6, 8, 9, 10, 11, 12, 13, 14]. The relation expresses a trade-off between
power fluctuations and entropy production. The entropy production imposes a
bound for the power fluctuations, suggesting that there may be fundamental
limitations on the precision of thermal machines. However, this is only a
bound for classical systems since quantum systems are able to circunvent the
TUR. The relevance of the TUR has been extended to be related with other
important concepts of nonequilibrium thermodynamics, including fluctuation
theorems [8, 9] and information theory [15, 16, 17].

Figure 1: (a) Scheme of the double quantum dot setup as a hybrid thermal
machine. Two quantum dots are capacitively coupled, and coupled to three
normal leads with $T_{r,h,c}$ and one superconducting lead. The heat current
from each lead to the system is denoted by $J_{r,h,c}$ and its sign criteria
is indicated. (b)Three working modes of this device (RP, ER, EP) are indicated
in the domain of $\mathbb{J}=(\dot{w},J_{h},J_{c})$, where
$\dot{w}=J_{r}+J_{h}+J_{c}$ means the production rate of work. For instance,
when a system is in RP mode (refrigerating and heat pumping),
$\dot{w}<0,J_{h}<0,J_{c}>0$. The EP mode is generating work and pumping heat,
$\dot{w}>0,J_{h}<0,J_{c}<0$. The heat current $J_{h,c}$ is positive when it
flows to the system.

Standard thermal machines consists of a working substance connected to two
heat baths, a cold reservoir and a hot one. Then, depending on the performed
task they can work as heat engines (E), refrigerators (R) or pumping machines
(P). Recently a novel characterization for systems connected to more than two
terminals, i.e., multiterminal devices, has been formulated. These are hybrid
thermal machines and have the advantage of performing more than one useful
task simultaneously. For these devices a generalized efficiency has been
enunciated [18]. Our goal is to extend the TUR for these kind of multiterminal
thermal machines and relate its efficiency (for the different working modes
listed below) with its current-current correlations. As previously stated, we
illustrate how these multidimensional TUR (MTUR) can be circumvented when the
working substance is a quantum system. For such case we show that MTUR are
violated and its cause is attributed to quantum coherence, to the non-local
character of quantum states in entangled states or the breakdown of LDB by
quantum effects [19]. Therefore, exploiting quantumness can be advantageous in
designing a precise quantum engine with low dissipation [20, 21, 22, 23] as
the generalized MTUR can be overcome. Definitely, the strong bound imposed by
the MTUR needs a revision for quantum systems where a careful analysis must be
done in view of coherence, interactions, or non-Markovian dynamics among
others [24, 25, 26, 27, 28, 29, 30].

We propose to deal with a minimal setup consisting in two nanoscale conductors
which are capacitively coupled, each of them connected to two baths (see Fig.
1). The setup is able to perform useful tasks, for instance, conducting an
electrical current, cooling a cold reservoir, or pumping heat when it is
electrically or thermal biased (or both). In this sense we are dealing with a
hybrid thermal machine [18]. Our study focus on the bounds obtained from the
thermodynamic uncertainty relations that are generalized to treat
multiterminal conductors. Departures from these bounds are illustrated with
our device in which the presence of coherent states and the breakdown of local
detailed balance (LDB) take place. We consider one of the contacts being a
superconductor, i.e., a reservoir of coherent Cooper pairs. Besides, each
conductor consists of a quantum dot, which is chosen due to its ability as
energy filters, a characteristic that potentially enhances its
thermoelectrical behavior.

_Hybrid thermal Machines: Efficiency and MTUR—_ Frequently, the efficiency of
a thermal engine is measured by the ratio of the output to the input. When
more than two terminals are attached to an engine, however, the physical
scenario becomes more involved. In general, multiterminal devices may exhibit
ill-defined efficiencies when they are characterized in a standard manner
[31]. The main cause is that the engine can produce multiple outputs
simultaneously. This occurs, for example, when work is extracted and at the
same time heat is pumped towards the hot reservoir [18] in three terminal
devices, the efficiency $\eta_{EP}$ is for a simultaneous engine-pumping mode.
Other operating modes correspond to refrigerator-pumping (RP) and engine-
refrigerator (ER) regimes. In the following, we will consider the situation in
which different outputs are produced concurrently, whereas inputs are heat
and/or consuming work. Quantifying and comparing the usefulness of different
outputs in response to the provided inputs has been discussed recently [32].
In there, the efficiency for a particular operating mode (consisting in a set
of inputs and outputs) is given as

$\eta^{\text{(mode)}}=\frac{\sum^{+}_{\alpha}\dot{w}^{\alpha}+\sum_{j}^{+}J_{j}\left(\frac{T_{\rm
r}}{T_{j}}-1\right)}{-\sum_{\alpha}^{-}\dot{w}^{\alpha}-\sum_{j}^{-}J_{j}\left(\frac{T_{\rm
r}}{T_{j}}-1\right)}\leq 1.$ (1)

where modes of hybrid machine $m\in\\{P,R,E,ER,EP,RP\\}$. Here,
$\sum^{\pm}_{j}J_{j}=\sum_{j}(J_{j}\pm|J_{j}|)/2$ are the sums over the
positive and negative heat currents $J_{j}$ and $\dot{w}$ is the production
rate of work performed by the system 111 The sign convention is adopted such
that a current is positive when it flows to the system. For instance, when a
system is in the RP mode, it pumps heat to a hot reservoir, $J_{h}<0$, and it
refrigerates a cold reservoir, $J_{c}>0$, by the work provided from the
environment, $\dot{w}<0$. In this case, the denominator of the efficiency
expression is $-\sum_{\alpha}^{-}\dot{w}^{\alpha}$, and the numerator is
$J_{j}\left(\frac{T_{\rm r}}{T_{j}}-1\right)$. That is, currents doing useful
works are put in the numerator, while a current providing a resource is put in
the denominator. . The choice of the reference temperature $T_{\rm r}$ depends
on the multiterminal configuration and the input/output tasks. Let us focus on
the case of three terminals. This scenario applies to our setup since although
the two quantum dot setup in Fig.1(a) contains the four leads, there is no
particle and energy flow to the superconducting (SC) lead in the large pairing
gap limit. The SC lead is analytically integrated in our treatment. (extension
to n-baths will be given elsewhere). The MTUR reads ($k_{B}=1$) [34]

$\mathbb{J}_{0}^{T}D^{-1}\mathbb{J}_{0}\leq\frac{1}{2}\sigma,$ (2)

where vector $\mathbb{J}_{0}=(J_{r},J_{h},J_{c})^{T}$ is a three component
vector with heat current via reference lead, hot lead, and cold lead,
respectively. The matrix element $D_{ij}$ contains the physical meaning of
thermal conductance between lead $i$ and $j$ according to the fluctuation-
dissipation theorem. The implication of the relation is that the precision of
measurement ($\mathbb{J}_{0}^{T}D^{-1}\mathbb{J}_{0}\equiv 1/\Delta P$) is
upper bounded by the entropy production $\sigma$. That is, the uncertainty of
measurement precision $\Delta P$ can be enhanced at the expense of entropy
production, $\sigma\,\Delta P\geq 2$.

To make a connection with the thermal efficiency, which is in terms of heat
current and the projection rate of work, we perform a transformation to
$\mathbb{J}=(J_{w},J_{h},J_{c})^{T}={\mathcal{A}}\mathbb{J}_{0}$ using
$J_{w}=\dot{w}=J_{r}+J_{h}+J_{c}=-\sum_{j=r,h,c}\mu_{j}I_{j}$, where
$\mathcal{A}$ is a 3$\times$3 matrix 222Explicitly,  $\displaystyle\mathbb{J}$
$\displaystyle=\left(\begin{array}[]{ccc}1&1&1\\\ 0&1&0\\\
0&0&1\end{array}\right)\left(\begin{array}[]{c}J_{r}\\\ J_{h}\\\
J_{c}\end{array}\right)=\mathcal{A}\mathbb{J}_{0}.$ (9) . The entropy
production rate is conveniently expressed by the current vector:

$\displaystyle\sigma$
$\displaystyle=\sum_{j=r,h,c}\frac{\dot{Q}_{j}}{T_{j}}=\frac{1}{T_{r}}\left(\mathbb{t_{w}}\cdot\mathbb{J}-\mathbb{t_{u}}\cdot\mathbb{J}\right),$
(10)

where $\mathbb{t}_{w}$ and $\mathbb{t}_{u}$ are vectors projecting the current
to the one generating and reducing the entropy, respectively. The former
current $J_{w}=\hat{\mathbb{t}}_{w}\cdot\mathbb{J}$ consumes a resource to
generate useful current $J_{u}=\hat{\mathbb{t}}_{u}\cdot\mathbb{J}$. For
instance, when a system is in the RP mode,
$\mathbb{t}_{w}^{\text{RP}}=(-1,0,0)$ and
$\mathbb{t}_{u}^{\text{RP}}=(0,-1+\frac{T_{r}}{T_{h}},\frac{T_{r}}{T_{c}}-1)$.
The machine receives work from the environment and then use the resource to
pump heat and to refrigerate the cold reservoir. The efficiency is expressed
as follows:

$\displaystyle\eta^{RP}=\frac{\mathbb{t}_{u}\cdot\mathbb{J}}{\mathbb{t}_{w}\cdot\mathbb{J}}=\eta^{R}+\eta^{P},$
(11)

where the hybrid efficiency is then divided into the two separate single
efficiencies. The other hybrid working modes can be similarly expressed by
choosing different $\mathbb{t}_{u,w}$ vectors while maintaining the same
vectorial expression for entropy production $\sigma$ and the efficiency
$\eta^{m}$ shown above 333For instance,
$\displaystyle\mathbb{t}_{w}^{\text{EP}}=(0,0,-\frac{T_{r}}{T_{c}}+1),\,\,\,\,\mathbb{t}_{u}^{\text{EP}}=(1,-1+\frac{T_{r}}{T_{h}},0),$
(12)
$\displaystyle\mathbb{t}_{w}^{\text{ER}}=(0,1-\frac{T_{r}}{T_{h}},0),\,\,\,\,\mathbb{t}_{u}^{\text{ER}}=(1,0,\frac{T_{r}}{T_{c}}-1).$
(13) where the ER mode is generating work and refrigerating, and the EP mode
is generating work and heat pumping.. We introduce one more vector
$\mathbb{t}^{\text{RP}}_{\perp}=(0,-\frac{T_{r}}{T_{c}}+1,-1+\frac{T_{r}}{T_{h}})$
to be perpendicular to both $\mathbb{t}_{w,u}$. The expression of the entropy
production and the efficiency motivate us to write the current vector in the
basis of $\mathbb{t}_{u,w,\perp}$.
$\mathbb{J}_{t}=(J_{w},J_{u},J_{\perp})^{T}=\mathcal{B}\mathbb{J},$ where
$\mathcal{B}$ is a 3$\times$3 matrix in terms of temperature $T_{r,h,c}$ only
444 For the RP operating mode of a hybrid machine,
$J_{u}=\left[J_{h}(-1+\frac{T_{r}}{T_{h}})+J_{c}(\frac{T_{r}}{T_{c}}-1)\right]/|\mathbb{t}_{u}|$
and
$J_{\perp}=\left[-J_{h}(\frac{T_{r}}{T_{h}}-1)+J_{c}(-1+\frac{T_{r}}{T_{c}})\right]/|\mathbb{t}_{\perp}|$,
where $|\mathbb{t}_{u}|=|\mathbb{t}_{\perp}|$. Hence,
$\displaystyle\mathcal{B}=\left(\begin{array}[]{ccc}1&0&0\\\
0&\cos\alpha&\sin\alpha\\\ 0&-\sin\alpha&\cos\alpha\end{array}\right),$ (17)
where $\cos\alpha=(-1+\frac{T_{r}}{T_{h}})/|\mathbb{t}_{u}|$ and
$\sin\alpha=(\frac{T_{r}}{T_{c}}-1)/|\mathbb{t}_{u}|$. . The MTUR is then
expressed in terms of new current vector:

$\mathbb{J}_{t}^{T}\tilde{D}^{-1}\mathbb{J}_{t}\leq\frac{1}{T_{r}}(\mathbb{t}_{u}\cdot\mathbb{J})(1-\eta^{-1}),$
(18)

where $\tilde{D}=[\mathcal{B}\mathcal{A}D\mathcal{A}^{T}\mathcal{B}^{T}]$. The
expression is utilized to obtain the upper bound of current of interest, for
example, ‘useful’ current, $J_{u}$, is upper bounded by

$\displaystyle
J_{u}\leq\frac{1}{T_{r}}|\mathbb{t}_{u}|(1-\eta^{-1})\left[\sum_{i,j=u,w,\perp}(\tilde{D}^{-1})_{ij}\frac{J_{i}J_{j}}{J_{u}^{2}}\right]^{-1}.$
(19)

where the ratio of currents $J_{i}/J_{u}$ is fixed once efficiency
$\eta^{R,P}$ are given 555In the RP mode,
$\frac{J_{w}}{J_{u}}=-(\eta^{RP})^{-1}$, and
$\frac{J_{\perp}}{J_{u}}=-\frac{1}{\eta^{RP}}\left[\eta^{P}r_{t}-\eta^{R}r_{t}^{-1}\right]$,
where $r_{t}=\frac{T_{r}/T_{c}-1}{-1+T_{r}/T_{h}}$. Hence, the upper bound in
Eq. (19) is determined by thermal conductance tensor, temperatures and
efficiency.  . The upper bound of the thermal efficiency is given by the
Clausius relation. For a given thermal efficiency, the MTUR provides the
quantitative upper bound of useful current in terms of temperatures and
thermal conductance tensor. That is practically relevant information when we
need the generation of a certain amount of current regardless of an optimal
thermal efficiency. For instance, at the optimal efficiency, it is well known
that the amount of useful current that can be generated is zero. When a system
contains quantum elements that promote the quantum coherence in particle and
energy transport, the upper bound of the MTUR can be exceeded, and as a result
one can harness resource from a hybrid thermal machine beyond the classical
upper bound. In the following we explicitly show in our double quantum dot
setup the situation where the MTUR relation is violated as a result of the LDB
breaking by Cooper pair transport.

Figure 2:  The hybrid thermal machine configured to show six working modes as
indicated in (a) by colors in the domain of
$(V_{\text{bias}},\varepsilon_{b})$. $\mu_{r}=-0.04$, $T_{h}=1.4$,
$T_{c}=0.6$, $\varepsilon_{t}=4$, $U_{t}=1$, $U_{tb}=2$. $\Gamma_{r}=0.05$,
$\Gamma_{S}=0.2$. The energies are indicated in unit of reference temperature,
$T_{r}=1$. (b) Shows the efficiency of corresponding hybrid modes. (c)
Illustrates the value of the MTUR, $\sigma\Delta P$. The regime with values
less than 2 shows the violation. (d) Shows the region which shows the highest
efficiency and the violation of MTUR simultaneously at
$\varepsilon_{r}/T_{r}=-3$.

_A superconducting hybrid machine—_ To illustrate the validity and the
departure of the MTUR due to quantum effects we consider parallel double
quantum dots labelled by ”t” (top) and ”b” (bottom). See Fig. 1(a) for the
scheme of the setup. The top dot is connected to normal and superconducting
leads and the intradot Coulomb interaction strength in this dot is denoted by
$U_{t}$. The bottom dot is attached to two normal leads with temperature
$T_{h,c}$. The two dots are coupled capacitively and the interaction strength
is $U_{tb}$. In the limit of large gap (i.e., $\Delta\to\infty$), the
superconducting lead can be easily integrated out [39]. Under this condition,
the Hamiltonian then reads ${\cal{H}}={\cal{H}}_{\rm C}^{\rm
eff}+{\cal{H}}_{R}+{\cal{H}}_{T}$ with [40]

$\displaystyle{\cal{H}}_{\rm C}^{\rm eff}$ $\displaystyle=$
$\displaystyle\varepsilon_{b}d_{b}^{\dagger}d_{b}+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\varepsilon_{t\sigma}d_{t\sigma}^{\dagger}d_{t\sigma}+\Gamma_{S}\left(d_{t\uparrow}^{\dagger}d_{t\downarrow}^{\dagger}+\rm{H.c.}\right)$
(20)
$\displaystyle+U_{tb}d_{b}^{\dagger}d_{b}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}d_{t\sigma}^{\dagger}d_{t\sigma}+U_{t}d_{t\uparrow}^{\dagger}d_{t\uparrow}d_{t\downarrow}^{\dagger}d_{t\downarrow}.$

Here, ${\cal{H}}_{R}$ is the Hamiltonian for the uncoupled leads. The
tunnelling from lead to dot is described by the Hamiltonian ${\cal{H}}_{T}$.
The hybridization between lead and dot is given by $\Gamma_{l}$ with
$l=\\{r,h,c\\}$. The energy levels of the top and bottom dots are indicated by
$\varepsilon_{t\sigma}$ and $\varepsilon_{b}$, respectively. The operators
$d_{t\sigma}^{\dagger}$ ($d_{t\sigma}$) and $d_{b}^{\dagger}$ ($d_{b}$) stand
for the creation (annihilation) operators for electrons in the top and bottom
dots with $\sigma=\\{\uparrow,\downarrow\\}$. Hereafter, we deploy our
formalism using the master equation approach. With this aim we diagonalize the
Hamiltonian given by Eq. (36). To do so, we take $|n_{b},n_{t}\rangle$ as the
basis, where $n_{b}=\\{0,1\\}$ is the occupation for the spinless bottom dot
and $n_{t}=\\{0,\uparrow/\downarrow,2\\}$ the corresponding occupation for the
top dot. When we diagonalize Eq. (36) in this basis, the hybridization
$\Gamma_{S}$ between the top dot and the superconducting lead mixes coherently
states with even number of electrons. The eigenstates are then given by
$|n_{b},\uparrow/\downarrow\rangle$, and
$|n_{b},\pm\rangle=u_{n_{b},\pm}|n_{b},0\rangle+v_{n_{b},\pm}|n_{b},2\rangle$,
where
$u_{n_{b},\pm}=(\varepsilon_{n_{b}}\pm\sqrt{\varepsilon_{n_{b}}^{2}+\Gamma_{S}^{2}})/\mathcal{N}_{n_{b},\pm}$
and $v_{n_{b},\pm}=-\Gamma_{S}/\mathcal{N}_{n_{b},\pm}$ with
$\varepsilon_{n_{b}}=\varepsilon_{t}+U_{t}/2+n_{b}U_{tb}$. Here,
$\mathcal{N}_{n_{b},\pm}$ is the normalization factor. The eigen energies for
$|n_{b},\pm\rangle$ are
$E_{\ket{n_{b}\pm}}=n_{b}\varepsilon_{b}+\varepsilon_{n_{b}}\mp\sqrt{\varepsilon_{n_{b}}^{2}+\Gamma_{S}^{2}}$
whereas for $|n_{b},\sigma\rangle$, these are simply
$E_{\ket{n_{b}\sigma}}=\varepsilon_{t}+n_{b}(\varepsilon_{b}+U_{tb})$.

_Pauli master equation—_ We apply the master equation formalism to describe
the dynamics for the occupation probability $\rho_{\theta}$ of the eigenstate
$\ket{\theta}$. To this end, we write down the transition rate from the
eigenstate $|\theta\rangle$ to $|\zeta\rangle$ as through the $\ell$ barrier
across the $X=t,b$ quantum dot

$\gamma_{l|\zeta\rangle\leftarrow|\theta\rangle}^{e(h)}=\Gamma_{l}\left|\langle\zeta|\Delta^{e(h)}_{X}|\theta\rangle\right|^{2}f^{e(h)}_{l}(E_{\theta}-E_{\zeta}),$
(21)

where $\Delta^{e(h)}_{t}=\sum_{\sigma}d_{t\sigma}^{\dagger}(d_{t\sigma})$ and
for $\Delta^{e(h)}_{b}=d^{\dagger}(d)$. In the limit
$k_{B}T_{l}\gg\Gamma_{l}$, the dynamics of the system is governed by the
sequential tunneling events and we safely neglect higher-order tunneling
correlations. As we are interested in the calculation of the current
fluctuations we write down the generalize master equation equation considering
the counting fields that generate the complete full counting statistics
$d\rho(\chi)/dt=-{\cal{L}}\rho(\chi$) where ${\cal{L}}$ means Louvillian
matrix. In the stationary limit ${\cal{L}}\rho=0$, we obtain the full counting
statistics formalism described in the supplemental material [41].

Figure 3: The quantity from which the violation of the TUR can be read off is
plotted for individual lead. (a) $\sigma\langle\Delta I_{r}^{2}\rangle/\langle
I_{r}\rangle^{2}$, associated with the charge current from the reference lead.
(b,c) $\sigma\langle\Delta J_{h,c}^{2}\rangle/\langle J_{h,c}\rangle^{2}$,
heat current from the hot and cold lead, respectively. (d)
$\sigma\langle\Delta\dot{w}^{2}\rangle/\langle\dot{w}\rangle^{2}$, associated
with the production rate of work (d). The transition between thermal working
modes is reflected in the divergence of the TUR value plotted. The breaking of
the local detailed balance appears in the particle current from the reference
lead due to the proximity effect from superconducting lead.

_Breakdown of Local Detailed Balance—_ Local detailed balance [42] amounts to
requiring that the log-ratio of an individual transition to the probability of
its time-reversed transition equals the entropy flux:
$\gamma_{F}/\gamma_{B}=e^{-\sigma_{C}}$. In systems where LDB is broken
exhibits different entropy productions when they are calculated by following
the Clausius relation ($\sigma_{C}$) and the Shannon entropy expression
($\sigma_{S}$). Due to the breakdown of the LDB an information flow ($I_{F}$)
is established by doing $\sigma_{S}=\sigma_{C}+I_{F}$ from the Maxwell demon
point of view [43, 44]. In our device LDB is manifestly broken having two
consequences, namely (i) the appearance of drag currents as reported in Ref.
[45], and (ii) the violation of the thermodynamic uncertainty relation [46]
which is our main focus. Due to the Andreev processes, the forward and
backward transition processes are not connected by Boltzmann weights. The
ratio of the forward ($\ket{n_{b},\sigma}\rightarrow\ket{n_{b},\pm}$) to
backward ($\ket{n_{b},\pm}\rightarrow\ket{n_{b},\sigma}$) transition is (See
[41] for full details)

$\frac{\gamma_{B}}{\gamma_{F}}=e^{\frac{\Delta
E}{T_{r}}}\left[\frac{\cosh\theta_{\gamma}+e^{-\frac{\Delta
E}{T_{r}}}\cosh(\theta_{\gamma}-\theta_{\mu_{r}})}{\cosh(\theta_{\gamma}+\theta_{\mu_{r}})+e^{-\frac{\Delta
E}{T_{r}}}\cosh\theta_{\gamma}}\right],$ (22)

where $\Delta E=E_{\ket{n_{b}\pm}}-E_{\ket{n_{b}\sigma}}$,
$\cosh\theta_{\gamma}=\frac{1}{2}(\left|\frac{v_{n_{b},\pm}}{u_{n_{b},\pm}}\right|^{2}+\left|\frac{u_{n_{b},\pm}}{v_{n_{b},\pm}}\right|^{2})$
and
$\cosh\theta_{\mu_{r}}=\frac{1}{2}(e^{\frac{\mu_{r}}{T_{r}}}+e^{-\frac{\mu_{r}}{T_{r}}})$.
$\mu_{r}$ is the chemical potential of reference lead. The ratio in Eq. (22)
depends on the the superconducting quasiparticle amplitudes
$u_{n_{b}\pm}/v_{n_{b}\pm}$ when $\mu_{r}\neq 0$, thus the LDB is broken. The
LDB is restored when the proximity effect in the top quantum dot is removed,
$\theta_{\gamma}\rightarrow\infty$. In brief, Andreev coherent processes are
responsible for the breakdown of the LDB under nonequilibrium conditions being
the main cause for the deviation of the TUR.

_Results: RP, EP, and ER modes—_ The versatile double quantum dot system shows
the all possible hybrid machine modes as indicated in Fig. 2(a) by tuning
$\varepsilon_{b}$ and $V_{\text{bias}}=\mu_{h}-\mu_{c}$. The hybrid machine
efficiency is achieved as high as $\eta\sim 0.85$ closely following the
boundary where the three machine modes (E, P and R) are present nearby each
other. For instance, the boundary of E and RP, that of R and EP, and that of P
and ER is the close neighbor of the high efficiency region shown in Fig. 2(b).
Note that the thermal efficiency expression is dependent of the modes, Eq.
(1). When the chemical potential of the reference lead $\mu_{r}$ is set
nonzero, the MTUR relation is violated at those boundaries of three machine
modes as shown in $\sigma\Delta P$ plotted in Fig. 2(c). This implies that it
is possible to achieve the optimally operating thermal machine generating
useful currents exceeding the upper bound limitation set by the classical
MTUR. The violation is associated with the quantum transport mainly involved
in the upper quantum within the proximity effect of superconducting lead. The
introduction of $\mu_{r}\neq 0$ triggers the quantum effect kicked in the
machine operation. The claim is supported in Fig. 3(a), where the TUR
associated with the charge current via the reference lead shows the violation,
$\sigma(\Delta P)_{r,part}\equiv\sigma\langle\Delta
I_{r}^{2}\rangle/I_{r}^{2}<2$.

The six hybrid machine modes are determined based on the sign of elements of
current vector $\mathbb{J}=(\dot{w},J_{h},J_{c})$. The mode of thermal machine
switches when a current element approaches to zero and then flips its sign. In
the course, the uncertainty of measurement precision $\Delta
P_{i}\equiv\langle\Delta J_{i}^{2}\rangle/J_{i}^{2},i\in\\{w,h,c\\}$ diverges.
Figure 3(b,c,d) shows $\sigma\Delta P_{i=h,c,w}$ associated with the current
of immediate interest, $J_{h,c,w}$. The phase boundary of the thermal machine
modes $P,R,E$ can be read off from the divergence of $\sigma\Delta
P_{i=h,c,w}$ which is consistent with the diagram of modes drawn in Fig. 2(a).
It is worth noting that the TUR of heat current associated with single lead,
$\sigma\Delta P_{i=r,h,c}$, does not show the violation, while the MTUR of
heat current drawn in Fig. 2(c) does. It implies that in a multiterminal
system the information contained in the correlation of cross terminal currents
is crucial to properly characterize the system.

Conclusion— The thermodynamic uncertainty relation has been formulated for
multiple-bath devices and it establishes a trade-off between generalized
current fluctuations and entropy production. A new class of machines, the so-
called hybrid machines, are devices that are able to perform several useful
tasks at once. Here, we have applied the uncertainty relation to deal with
multiterminal devices working as hybrid machines. From the MTUR a new bound is
established for the resourceful current in terms of the generalized efficiency
and the cross current correlations. We have shown that MTUR are violated in
quantum systems as a result of the LDB breaking. We have shown that depending
on the operating mode (RP, EP, or ER) a violation of the MTUR is accompanied
by a high efficiency being the optimal scenario for an hybrid thermal machine.
We provide a general guideline for the design of multiterminal hybrid machines
in the quantum regime that operate at the best task performance, meaning low
dissipation, small fluctuations around average values and high task efficiency
with a moderately high resourceful current.

Acknowledgements— R.L acknowledges the financial support through the grants
PID2020-117347GB-I00 and the grant from the María de Maeztu Program for Units
of Excellence No. MDM2017-0711 funded by MCIN/AEI/10.13039/501100011033’.
K.W.K. acknowledges financial support by Basic Science Research Program
through the National Research Foundation of Korea (NRF) funded by the Ministry
of Education (20211060) and Korea government(MSIT) (No.2020R1A5A1016518).

## References

  * Carnot [1978] S. Carnot,  _Réflexions sur la puissance motrice du feu_ , 26 (Vrin, 1978).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic uncertainty relation for biomolecular processes, Physical review letters 114, 158101 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Dissipation bounds all steady-state current fluctuations, Physical review letters 116, 120601 (2016).
  * Pietzonka  _et al._ [2016a] P. Pietzonka, A. C. Barato, and U. Seifert, Universal bounds on current fluctuations, Physical Review E 93, 052145 (2016a).
  * Pietzonka  _et al._ [2016b] P. Pietzonka, A. C. Barato, and U. Seifert, Universal bound on the efficiency of molecular motors, Journal of Statistical Mechanics: Theory and Experiment 2016, 124004 (2016b).
  * Pietzonka and Seifert [2018] P. Pietzonka and U. Seifert, Universal trade-off between power, efficiency, and constancy in steady-state heat engines, Physical review letters 120, 190602 (2018).
  * Gingrich and Horowitz [2017] T. R. Gingrich and J. M. Horowitz, Fundamental bounds on first passage time fluctuations for currents, Physical review letters 119, 170601 (2017).
  * Timpanaro  _et al._ [2019] A. M. Timpanaro, G. Guarnieri, J. Goold, and G. T. Landi, Thermodynamic uncertainty relations from exchange fluctuation theorems, Physical review letters 123, 090604 (2019).
  * Hasegawa and Van Vu [2019] Y. Hasegawa and T. Van Vu, Fluctuation theorem uncertainty relation, Physical review letters 123, 110602 (2019).
  * Pietzonka  _et al._ [2017] P. Pietzonka, F. Ritort, and U. Seifert, Finite-time generalization of the thermodynamic uncertainty relation, Physical Review E 96, 012101 (2017).
  * Horowitz and Gingrich [2017] J. M. Horowitz and T. R. Gingrich, Proof of the finite-time thermodynamic uncertainty relation for steady-state currents, Physical Review E 96, 020103 (2017).
  * Shiraishi [2017] N. Shiraishi, Finite-time thermodynamic uncertainty relation do not hold for discrete-time markov process, arXiv preprint arXiv:1706.00892 (2017).
  * Proesmans and Van den Broeck [2017] K. Proesmans and C. Van den Broeck, Discrete-time thermodynamic uncertainty relation, EPL (Europhysics Letters) 119, 20001 (2017).
  * Macieszczak  _et al._ [2018] K. Macieszczak, K. Brandner, and J. P. Garrahan, Unified thermodynamic uncertainty relations in linear response, Physical review letters 121, 130601 (2018).
  * Agarwalla and Segal [2018a] B. K. Agarwalla and D. Segal, Assessing the validity of the thermodynamic uncertainty relation in quantum systems, Physical Review B 98, 155438 (2018a).
  * Marsland and England [2017] R. Marsland and J. England, Limits of predictions in thermodynamic systems: a review, Reports on Progress in Physics 81, 016601 (2017).
  * Potts and Samuelsson [2019] P. P. Potts and P. Samuelsson, Thermodynamic uncertainty relations including measurement and feedback, Physical Review E 100, 052137 (2019).
  * Manzano  _et al._ [2020a] G. Manzano, R. Sánchez, R. Silva, G. Haack, J. B. Brask, N. Brunner, and P. P. Potts, Hybrid thermal machines: Generalized thermodynamic resources for multitasking, Phys. Rev. Research 2, 043302 (2020a).
  * López  _et al._ [2012] R. López, J. S. Lim, and D. Sánchez, Fluctuation relations for spintronics, Phys. Rev. Lett. 108, 246603 (2012).
  * Ptaszyński [2018] K. Ptaszyński, Coherence-enhanced constancy of a quantum thermoelectric generator, Physical Review B 98, 085425 (2018).
  * Agarwalla and Segal [2018b] B. K. Agarwalla and D. Segal, Assessing the validity of the thermodynamic uncertainty relation in quantum systems, Physical Review B 98, 155438 (2018b).
  * Brandner  _et al._ [2018a] K. Brandner, T. Hanazato, and K. Saito, Thermodynamic bounds on precision in ballistic multiterminal transport, Physical review letters 120, 090601 (2018a).
  * Liu and Segal [2019] J. Liu and D. Segal, Thermodynamic uncertainty relation in quantum thermoelectric junctions, Physical Review E 99, 062141 (2019).
  * Michał Horodecki [2013] J. O. Michał Horodecki, Fundamental limitations for quantum and nanoscale thermodynamics, Nature Communications 4, 2059 (2013).
  * Segal [2018] D. Segal, Current fluctuations in quantum absorption refrigerators, Physical Review E 97, 052145 (2018).
  * Brandner  _et al._ [2018b] K. Brandner, T. Hanazato, and K. Saito, Thermodynamic bounds on precision in ballistic multiterminal transport, Phys. Rev. Lett. 120, 090601 (2018b).
  * Ptaszyński [2018] K. Ptaszyński, Coherence-enhanced constancy of a quantum thermoelectric generator, Phys. Rev. B 98, 085425 (2018).
  * Saryal  _et al._ [2019] S. Saryal, H. M. Friedman, D. Segal, and B. K. Agarwalla, Thermodynamic uncertainty relation in thermal transport, Physical Review E 100, 042101 (2019).
  * Buffoni and Campisi [2020] L. Buffoni and M. Campisi, Thermodynamics of a quantum annealer, Quantum Science and Technology 5, 035013 (2020).
  * Benenti  _et al._ [2020] G. Benenti, G. Casati, and J. Wang, Power, efficiency, and fluctuations in steady-state heat engines, Physical Review E 102, 040103 (2020).
  * Hajiloo  _et al._ [2020] F. Hajiloo, R. Sánchez, R. S. Whitney, and J. Splettstoesser, Quantifying nonequilibrium thermodynamic operations in a multiterminal mesoscopic system, Physical Review B 102, 155405 (2020).
  * Manzano  _et al._ [2020b] G. Manzano, R. Sánchez, R. Silva, G. Haack, J. B. Brask, N. Brunner, and P. P. Potts, Hybrid thermal machines: Generalized thermodynamic resources for multitasking, Physical Review Research 2, 043302 (2020b).
  * Note [1] The sign convention is adopted such that a current is positive when it flows to the system. For instance, when a system is in the RP mode, it pumps heat to a hot reservoir, $J_{h}<0$, and it refrigerates a cold reservoir, $J_{c}>0$, by the work provided from the environment, $\dot{w}<0$. In this case, the denominator of the efficiency expression is $-\sum@\slimits@_{\alpha}^{-}\dot{w}^{\alpha}$, and the numerator is $J_{j}\left(\frac{T_{\rm r}}{T_{j}}-1\right)$. That is, currents doing useful works are put in the numerator, while a current providing a resource is put in the denominator.
  * Dechant [2018] A. Dechant, Multidimensional thermodynamic uncertainty relations, Journal of Physics A: Mathematical and Theoretical 52, 035001 (2018).
  * Note [2] Explicitly,  $\displaystyle\mathbf{J}$ $\displaystyle=\left(\begin{array}[]{ccc}1&1&1\\\ 0&1&0\\\ 0&0&1\end{array}\right)\left(\begin{array}[]{c}J_{r}\\\ J_{h}\\\ J_{c}\end{array}\right)=\mathcal{A}\mathbf{J}_{0}.$ (29) .
  * Note [3] For instance,  $\displaystyle\mathbf{t}_{w}^{\text{EP}}=(0,0,-\frac{T_{r}}{T_{c}}+1),\,\,\,\,\mathbf{t}_{u}^{\text{EP}}=(1,-1+\frac{T_{r}}{T_{h}},0),$ (30) $\displaystyle\mathbf{t}_{w}^{\text{ER}}=(0,1-\frac{T_{r}}{T_{h}},0),\,\,\,\,\mathbf{t}_{u}^{\text{ER}}=(1,0,\frac{T_{r}}{T_{c}}-1).$ (31) where the ER mode is generating work and refrigerating, and the EP mode is generating work and heat pumping.
  * Note [4] For the RP operating mode of a hybrid machine, $J_{u}=\left[J_{h}(-1+\frac{T_{r}}{T_{h}})+J_{c}(\frac{T_{r}}{T_{c}}-1)\right]/|\mathbf{t}_{u}|$ and $J_{\perp}=\left[-J_{h}(\frac{T_{r}}{T_{h}}-1)+J_{c}(-1+\frac{T_{r}}{T_{c}})\right]/|\mathbf{t}_{\perp}|$, where $|\mathbf{t}_{u}|=|\mathbf{t}_{\perp}|$. Hence,  $\displaystyle\mathcal{B}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&\mathop{cos}\nolimits\alpha&\mathop{sin}\nolimits\alpha\\\ 0&-\mathop{sin}\nolimits\alpha&\mathop{cos}\nolimits\alpha\end{array}\right),$ (35) where $\mathop{cos}\nolimits\alpha=(-1+\frac{T_{r}}{T_{h}})/|\mathbf{t}_{u}|$ and $\mathop{sin}\nolimits\alpha=(\frac{T_{r}}{T_{c}}-1)/|\mathbf{t}_{u}|$.
  * Note [5] In the RP mode, $\frac{J_{w}}{J_{u}}=-(\eta^{RP})^{-1}$, and $\frac{J_{\perp}}{J_{u}}=-\frac{1}{\eta^{RP}}\left[\eta^{P}r_{t}-\eta^{R}r_{t}^{-1}\right]$, where $r_{t}=\frac{T_{r}/T_{c}-1}{-1+T_{r}/T_{h}}$. Hence, the upper bound in Eq. (19) is determined by thermal conductance tensor, temperatures and efficiency.
  * Rozhkov and Arovas [1999] A. V. Rozhkov and D. P. Arovas, Josephson coupling through a magnetic impurity, Phys. Rev. Lett. 82, 2788 (1999).
  * Rozhkov and Arovas [2000] A. Rozhkov and D. P. Arovas, Interacting-impurity josephson junction: Variational wave functions and slave-boson mean-field theory, Physical Review B 62, 6687 (2000).
  * [41] Supplemental material, supplemental Material.
  * Katz  _et al._ [1983] S. Katz, J. L. Lebowitz, and H. Spohn, Phase transitions in stationary nonequilibrium states of model lattice systems, Phys. Rev. B 28, 1655 (1983).
  * Esposito and Schaller [2012] M. Esposito and G. Schaller, Stochastic thermodynamics for “maxwell demon” feedbacks, EPL (Europhysics Letters) 99, 30003 (2012).
  * Rosselló  _et al._ [2017] G. Rosselló, R. López, and G. Platero, Chiral maxwell demon in a quantum hall system with a localized impurity, Phys. Rev. B 96, 075305 (2017).
  * Tabatabaei  _et al._ [2020] S. M. Tabatabaei, D. Sánchez, A. L. Yeyati, and R. Sánchez, Andreev-coulomb drag in coupled quantum dots, Phys. Rev. Lett. 125, 247701 (2020).
  * Maes [2021] C. Maes, Local detailed balance, SciPost Phys. Lect. Notes , 32 (2021).
  * Esposito and Mukamel [2006] M. Esposito and S. Mukamel, Fluctuation theorems for quantum master equations, Physical Review E 73, 046129 (2006).
  * Shiraishi and Saito [2019] N. Shiraishi and K. Saito, Fundamental relation between entropy production and heat current, Journal of Statistical Physics 174, 433 (2019).
  * Maisi  _et al._ [2014] V. F. Maisi, D. Kambly, C. Flindt, and J. P. Pekola, Full counting statistics of andreev tunneling, Phys. Rev. Lett. 112, 036801 (2014).
  * Flindt  _et al._ [2008] C. Flindt, T. c. v. Novotný, A. Braggio, M. Sassetti, and A.-P. Jauho, Counting statistics of non-markovian quantum stochastic processes, Phys. Rev. Lett. 100, 150601 (2008).
  * Levitov  _et al._ [1996] L. S. Levitov, H. Lee, and G. B. Lesovik, Electron counting statistics and coherent states of electric current, Journal of Mathematical Physics 37, 4845 (1996).
  * Lesovik and Chtchelkatchev [2003] G. Lesovik and N. Chtchelkatchev, Quantum and classical binomial distributions for the charge transmitted through coherent, JEPT Lett. 77, 393 (2003).
  * Drazin [1958] M. P. Drazin, Pseudo-inverses in associative rings and semigroups, The American Mathematical Monthly 65, 506 (1958).

##  I Supplemental Materials

###### Contents

  1. I Supplemental Materials
  2. II A double quantum dot system
  3. III Entropy production: Shannon and Clausius
  4. IV The local detailed balance
  5. V Hybrid thermal machine
  6. VI Full Counting Statistics
  7. VII Quantum master equation
    1. VII.1 A single quantum dot
    2. VII.2 A double quantum dot

##  II A double quantum dot system

In this section we provide the detailed expressions of eigenenergies and
corresponding eigenmodes of the double quantum dot system. We allow the double
occupancy in the top quantum dot, and single occupancy in the bottom quantum
dot. There is an onsite interaction in the top quantum dot ($U_{t}$), and an
intra-dot interaction ($U_{tb}$) is present by a capacitive coupling. The
bottom dot is coupled to a hot reservoir and cold reservoir, and their
chemical potential difference ($\mu_{h}=V_{\rm bias}/2$ and $\mu_{c}=-V_{\rm
bias}/2$) makes a electrical bias across the bottom quantum dot. The top
quantum dot is coupled to one metallic reference lead with temperature $T_{r}$
and one superconducting lead. The latter provides the superconducting
proximity effect. We consider the regime where the system-lead coupling is
small, and they will be treated perturbatively. First, the Hamiltonian of the
double quantum dot system is:

$\displaystyle{\cal{H}}_{\rm DQD}^{\rm eff}$ $\displaystyle=$
$\displaystyle\varepsilon_{b}d_{b}^{\dagger}d_{b}+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\varepsilon_{t}d_{t\sigma}^{\dagger}d_{t\sigma}+\Gamma_{S}\left(d_{t\uparrow}^{\dagger}d_{t\downarrow}^{\dagger}+\rm{H.c.}\right)$
(36)
$\displaystyle+U_{tb}d_{b}^{\dagger}d_{b}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}d_{t\sigma}^{\dagger}d_{t\sigma}+U_{t}d_{t\uparrow}^{\dagger}d_{t\uparrow}d_{t\downarrow}^{\dagger}d_{t\downarrow}.$

where $\varepsilon_{t,b}$ is the chemical potential of top and bottom quantum
dot, respectively. $\Gamma_{S}$ is the coupling strength between the top
quantum dot and the superconducting lead. To diagonalize ${\cal{H}}_{\rm
DQD}^{\rm eff}$, we consider the basis,
$|n_{b},n_{t\uparrow},n_{t\downarrow}\rangle=\left(d_{b}^{\dagger}\right)^{n_{b}}\left(d_{t\uparrow}^{\dagger}\right)^{n_{t\uparrow}}\left(d_{t\downarrow}^{\dagger}\right)^{n_{t\downarrow}}|0\rangle$,
where $n_{b}=0,1$ and $n_{t\sigma}=0,1$. Thus, there are total eight states in
the double quantum dot system. Specifically,

Figure 4: The entropy production rate according to the Shannon expression
Eq.(55) (left panel), $\sigma_{S}$. The difference between two expressions of
entropy production rate, $\sigma_{C}-\sigma_{S}$, is plotted (center panel).
When the local detailed balance is present, $\mu_{r}=0$, the difference
$\sigma_{C}-\sigma_{S}$ is plotted (right panel), showing that
$\sigma_{S}=\sigma_{C}$ within a numerical error.

$\displaystyle|0,\uparrow\rangle=|0,1,0\rangle,\quad\text{with}\quad
E_{\ket{0,\uparrow}}=\varepsilon_{t}$ (37)
$\displaystyle|0,\downarrow\rangle=|0,0,1\rangle,\quad\text{with}\quad
E_{\ket{0,\downarrow}}=\varepsilon_{t}$ (38)
$\displaystyle|1,\uparrow\rangle=|1,1,0\rangle,\quad\text{with}\quad
E_{\ket{1,\uparrow}}=\varepsilon_{b}+\varepsilon_{t}+U_{tb}$ (39)
$\displaystyle|1,\downarrow\rangle=|1,0,1\rangle,\quad\text{with}\quad
E_{\ket{1,\downarrow}}=\varepsilon_{b}+\varepsilon_{t}+U_{tb}$ (40)
$\displaystyle|0,+\rangle$
$\displaystyle=u_{0,+}|0,0,0\rangle+v_{0,+}|0,1,1\rangle,$
$\displaystyle\quad\quad\quad\text{with}\quad
E_{\ket{0,+}}=\varepsilon_{0}-\sqrt{\varepsilon_{0}^{2}+\Gamma_{S}^{2}},$ (41)
$\displaystyle|0,-\rangle$
$\displaystyle=u_{0,-}|0,0,0\rangle+v_{0,-}|0,1,1\rangle,$
$\displaystyle\quad\quad\quad\text{with}\quad
E_{\ket{0,-}}=\varepsilon_{0}+\sqrt{\varepsilon_{0}^{2}+\Gamma_{S}^{2}},$ (42)
$\displaystyle|1,+\rangle$
$\displaystyle=u_{1,+}|1,0,0\rangle+v_{1,+}|1,1,1\rangle,$
$\displaystyle\quad\quad\quad\text{with}\quad
E_{\ket{1,+}}=\varepsilon_{b}+\varepsilon_{1}-\sqrt{\varepsilon_{1}^{2}+\Gamma_{S}^{2}},$
(43) $\displaystyle|1,-\rangle$
$\displaystyle=u_{1,-}|1,0,0\rangle+v_{1,-}|1,1,1\rangle,$
$\displaystyle\quad\quad\quad\text{with}\quad
E_{\ket{1,-}}=\varepsilon_{b}+\varepsilon_{1}+\sqrt{\varepsilon_{1}^{2}+\Gamma_{S}^{2}},$
(44)

where we introduce the notation used in the main text. The coefficients of the
superconducting states are:

$u_{n_{b},\pm}=\frac{\varepsilon_{n_{b}}\pm\sqrt{\varepsilon_{n_{b}}^{2}+\Gamma_{S}^{2}}}{{\cal{N}}_{n_{b},\pm}},\quad
v_{n_{b},\pm}=-\frac{\Gamma_{S}}{{\cal{N}}_{n_{b},\pm}},$ (45)

with $u_{n_{b},\pm}^{2}+v_{n_{b},\pm}^{2}=1$. And,

$\varepsilon_{n_{b}}=\varepsilon_{t}+U_{p}/2+n_{b}U_{tb}\quad\text{for
{\hbox{n_b=0,1}}}\quad$ (46)

The normalization factor is given by

${\cal{N}}_{n_{b},\pm}=\sqrt{\left(\varepsilon_{n_{b}}\pm\sqrt{\varepsilon_{n_{b}}^{2}+\Gamma_{S}^{2}}\right)^{2}+\Gamma_{S}^{2}}$
(47)

For later use, we specify the energy difference between states by the
transition made by a coupled lead when $\Gamma_{S}\ll
U_{t},U_{tb},\varepsilon_{t}$, which is the working regime of our study:
$\Gamma_{S}=0.2$, $U_{t}=1$, $U_{tb}=2$, $\varepsilon_{t}=4$, and
$\varepsilon_{b}\in[-5,1]$ in the main text. The transitions between the eight
states are induced by the leads (reference, hot, and cold reservoirs). When
the transition involves an energy difference, the leads provide or absorb the
extra energy. In particular, the transitions induced by hot or cold lead are
the following (they make change in the particle number sitting in the bottom
quantum dot, $n_{b}$):

$\displaystyle E_{1,-}-E_{0,+}$
$\displaystyle\simeq\varepsilon_{b}+2\varepsilon_{t}+U_{t}+2U_{tb},$ (48)
$\displaystyle E_{1,-}-E_{0,-}$ $\displaystyle\simeq\varepsilon_{b}+2U_{tb},$
(49) $\displaystyle E_{1,\uparrow}-E_{0,\uparrow}$
$\displaystyle\simeq\varepsilon_{b}+U_{tb},$ (50) $\displaystyle
E_{1,\downarrow}-E_{0,\downarrow}$
$\displaystyle\simeq\varepsilon_{b}+U_{tb},$ (51) $\displaystyle
E_{1,+}-E_{0,+}$ $\displaystyle\simeq\varepsilon_{b},$ (52) $\displaystyle
E_{1,+}-E_{0,-}$ $\displaystyle\simeq\varepsilon_{b}-2\varepsilon_{t}-U_{t},$
(53)

where the energy difference is written in such a way that from the top to
bottom the energy difference decreases. They all share the energy
$\varepsilon_{b}$, in addition to the onsite and inter-dot interaction energy
$U_{t,tb}$. Later, we introduce the energy dependent system-(hot,cold) lead
coupling strength to operate the hybrid thermal machine, see Sec.V.

##  III Entropy production: Shannon and Clausius

The entropy is the number of accessible states in reservoirs, which is in our
case the three leads. When heat flows into a reservoir, the particles in the
reservoir can access more number of states and as a result the entropy
increases. The entropy production therefore can be expressed in terms of heat
currents, and this is the Clausius entropy $\sigma_{C}$:

$\displaystyle\sigma_{C}=\sum_{l=r,h,c}-\frac{J_{l}}{T_{l}},$ (54)

where the (-) sign is because we adopt the sign convention of current in such
way that it is positive when it flows to the system.

On the other and, the Shannon entropy production is obtained from the
transition rates between states of system, which are specified in Sec.II, and
the probability density [47]:

$\displaystyle\sigma_{S}=-\sum_{m,n}\gamma_{mn}\rho_{n}\ln\left(\frac{\gamma_{mn}\rho_{n}}{\gamma_{nm}\rho_{m}}\right).$
(55)

The two expressions of the entropy production in identical in the presence of
the local detailed balance (for example, see Sec.2.2.4 of Ref. [48]). In our
double quantum dot system, the local detailed balance can be broken by
introducing nonzero chemical potential $\mu_{r}$ at the reference frame, which
is coupled to the top quantum dot (See Sec.IV). In such a case the Clausius
($\sigma_{C}\equiv\sigma_{Clausius}$) and Shannon entropies
($\sigma_{S}\equiv\sigma_{Shannon}$) differ as shown in Fig. 4 for our setup.

Figure 5: Thermal efficiency (left) and MTUR value $\sigma\Delta P$ (right)
are plotted at chemical potential of the reference reservoir
$\mu_{r}/T_{r}=0,-0.01,-0.02,-0.03$ (indicated in the legend).

##  IV The local detailed balance

The local detailed balance (LDB) is broken in the transition which involves
the change of particle number in the top quantum dot when $\mu_{r}\neq 0$. The
transition is induced by the reference lead. Below, we consider the transition
from $\ket{n_{b}=0,\sigma\in\\{\uparrow,\downarrow\\}}$ to $\ket{n_{b}=0,+}$.
The similar argument follows for the transition
$\ket{n_{b}=1,\sigma\in\\{\uparrow,\downarrow\\}}$ to $\ket{n_{b}=1,\pm}$. The
forward transition rate is (the transition rate below is shown in Eq. (123)):

$\displaystyle\gamma^{r}_{\ket{0,+}\leftarrow\ket{0,\sigma}}$
$\displaystyle=\gamma^{(e)}_{\ket{0,+}\leftarrow\ket{0,\sigma}}+\gamma^{(h)}_{\ket{0,+}\leftarrow\ket{0,\sigma}},$
$\displaystyle=\Gamma_{r}\left(|v_{0+}|^{2}f^{e}(\Delta
E)+|u_{0+}|^{2}f^{h}(-\Delta E)\right),$
$\displaystyle=\Gamma_{r}\left[\frac{|v_{0+}|^{2}}{e^{(\Delta
E-\mu_{r})/T_{r}}+1}+\frac{|u_{0+}|^{2}}{e^{(\Delta
E+\mu_{r})/T_{r}}+1}\right],$ (56)

where $\Delta E=E_{\ket{0,+}}-E_{\ket{0,\sigma}}$, $f^{e}(\Delta
E)=(e^{(\Delta E-\mu_{r})/T_{r}}+1)^{-1}$, and $f^{h}(-\Delta
E)=1-(e^{(-\Delta E-\mu_{r})/T_{r}}+1)^{-1}$. (Energies are written in unit of
temperature of reference reservoir, $T_{r}$). There are two scattering
processes: one is to add electron to the top quantum dot that comes with
probability $|v_{0+}|^{2}$ and distribution $f^{e}(\Delta E)$, and the second
term is to remove an electron from the top quantum dot that comes with
probability $|u_{0+}|^{2}$ and distribution $f^{h}(-\Delta E)=1-f^{e}(-\Delta
E)$. For the second process, an electron is going into the reference lead with
the energy $\Delta E$ below $\mu_{r}$. The forward transition rate is
similarly written (the transition rate below is shown in Eq. (126)):

$\displaystyle\gamma^{r}_{\ket{0,\sigma}\leftarrow\ket{0,+}}$
$\displaystyle=\gamma^{(e)}_{\ket{0,\sigma}\leftarrow\ket{0,+}}+\gamma^{(h)}_{\ket{0,\sigma}\leftarrow\ket{0,+}},$
$\displaystyle=\Gamma_{r}\left(|v_{0+}|^{2}f^{e}(-\Delta
E)+|u_{0+}|^{2}f^{h}(\Delta E)\right),$
$\displaystyle=\Gamma_{r}\left[\frac{|v_{0+}|^{2}}{e^{(-\Delta
E-\mu_{r})/T_{r}}+1}+\frac{|u_{0+}|^{2}}{e^{(-\Delta
E+\mu_{r})/T_{r}}+1}\right],$ (57)

The ratio between the two transition rate is generally dependent of
$|v_{0+}|^{2}$ and $|u_{0+}|^{2}$ for $\mu_{r}\neq 0$, therefore the local
detailed balance is broken:

$\displaystyle\frac{\gamma_{B}}{\gamma_{F}}$
$\displaystyle=\frac{\gamma^{r}_{\ket{0,\sigma}\leftarrow\ket{0,+}}}{\gamma^{r}_{\ket{0,+}\leftarrow\ket{0,\sigma}}},$
$\displaystyle=\frac{|v_{0+}|^{2}f^{e}(-\Delta E)+|u_{0+}|^{2}f^{h}(\Delta
E)}{|v_{0+}|^{2}f^{e}(\Delta E)+|u_{0+}|^{2}f^{h}(-\Delta E)},$
$\displaystyle=e^{\frac{\Delta
E}{T_{r}}}\left[\frac{\cosh\theta_{\gamma}+e^{-\frac{\Delta
E}{T_{r}}}\cosh(\theta_{\gamma}-\theta_{\mu_{r}})}{\cosh(\theta_{\gamma}+\theta_{\mu_{r}})+e^{-\frac{\Delta
E}{T_{r}}}\cosh\theta_{\gamma}}\right],$ (58)

where
$\cosh\theta_{\gamma}=\frac{1}{2}(\left|\frac{v_{n_{b},\pm}}{u_{n_{b},\pm}}\right|^{2}+\left|\frac{u_{n_{b},\pm}}{v_{n_{b},\pm}}\right|^{2})$
and
$\cosh\theta_{\mu_{r}}=\frac{1}{2}(e^{\frac{\mu_{r}}{T_{r}}}+e^{-\frac{\mu_{r}}{T_{r}}})$.
If $\mu_{r}=0$, the local detailed balance is preserved. Since
$\theta_{\mu_{r}}=0$,

$\displaystyle\frac{\gamma_{B}}{\gamma_{F}}=\left.\frac{\gamma_{\ket{0,\sigma}\leftarrow\ket{0,+}}}{\gamma_{{\ket{0,+}\leftarrow\ket{0,\sigma}}}}\right|_{\mu_{r}=0}=e^{\Delta
E/T_{r}}.$ (59)

This implies that for $\Delta E>0$ the transition rate of the forward
transition is exponentially small in $\Delta E$ than that of the backward
transition. Due to the breakdown of the local detailed balance when
$\mu_{r}\neq 0$ the multidimensional thermodynamic uncertainty relation (MTUR)
is not satisfied as shown in Fig. 5 where the efficiency (left figure) and the
MTUR (right figure) are shown. For $\mu_{r}\in[0.01,0.05]$ the MTUR are
violated whereas the efficiency remains almost untouchable (see left figure in
Fig. 5).

Note that the role of superconducting lead is crucial for the breakdown of
LDB. Due to the proximity effect, the two states in Eq.(II-44) with even
number of particle difference are coupled. A situation without
superconductivity can be simulated by setting either $v_{0+}=0$ or $u_{0+}=0$.
In this case, $\theta_{\gamma}\rightarrow\infty$ in Eq.(IV) and the LDB is
recovered regardless of the chemical potential in the reference lead.

Figure 6: Hybrid machine mode (first row), thermal efficiency (second row),
and $\sigma\Delta P$ (third row) for different configuration of of energy-
dependent lead-system coupling. (Left column): $\gamma_{0,1,2}=\Gamma_{r}$
corresponding to the energy independent case. (Center column):
$\gamma_{0,2}=\Gamma_{r}$ and $\gamma_{1}=0$. (Right column):
$\gamma_{0}=\Gamma_{r}$ and $\gamma_{1,2}=0$.

##  V Hybrid thermal machine

The double quantum dot system is coupled to the four lead (reference,
superconducting, hot, and cold reservoir). Only three of them (except the
superconducting lead) can exchange particle and energy current with the system
because the superconducting gap is assumed to be infinite. Thus, our system is
effectively considered as a hybrid thermal machine that is capable of one or
two among three different machine modes, engine, heat pumping, and
refrigerating. When the system-lead coupling strength is set energy-
independent value, $\Gamma_{r,h,c}=0.05$, we find the working mode, thermal
efficiency, and $\sigma\Delta P$ as shown in the left column of Fig. 6. The
hybrid machine modes ER and EP are not present. To access the modes, it turns
out [32] that it is necessary to introduce energy-dependent system-lead
coupling for hot and cold reservoir.

$\displaystyle\Gamma_{h}(\Delta\tilde{E}\geq 3)=\gamma_{2},\quad$
$\displaystyle\Gamma_{c}(\Delta\tilde{E}\geq 3)=\gamma_{0},$ (60)
$\displaystyle\Gamma_{h}(3>\Delta\tilde{E}\geq 1)=\gamma_{1},\quad$
$\displaystyle\Gamma_{c}(3>\Delta\tilde{E}\geq 1)=\gamma_{1},$ (61)
$\displaystyle\Gamma_{h}(1>\Delta\tilde{E})=\gamma_{0},\quad$
$\displaystyle\Gamma_{c}(1>\Delta\tilde{E})=\gamma_{2},$ (62)

where ${\Delta\tilde{E}}=E_{\ket{f}}-E_{\ket{i}}-\varepsilon_{b}$ (see
Eq.(48-53) for the relevant transitions). For instance, we obtain the results
shown in the main text by setting $\gamma_{2}=\Gamma_{r},\gamma_{0,1}=0$. It
implies that the hot lead only exchange particles with energy
$\Delta\tilde{E}\geq 3$, while the cold lead only exchange particles with
energy $\Delta\tilde{E}<0$. Hence, when there is a particle flows from the hot
to cold lead, the energy current is correlated with it. In the center column
of Fig. 6, we set $\gamma_{0,2}=\Gamma_{r}$ and $\gamma_{1}=0$. In the right
column, we set $\gamma_{0}=\Gamma_{r}$ and $\gamma_{1,2}=0$, which is the
opposite situation from the condition in the main text: energetic particles
$\Delta\tilde{E}>2$ from the cold lead can tunnel to the system, while only
particles with $\Delta\tilde{E}<0$ from the hot lead can tunnel to the system.

##  VI Full Counting Statistics

The FCS technique is considered to derive the second cumulant for the heat and
charge transport and its correlations. For such purpose we follow Ref. [49,
50]. We firstly solve the set of master equations and obtain the cumulant
generating function [51, 52]. We first, introduce the eigenvectors $a_{\beta}$
and eigenvalues $\lambda_{\beta}$ of the matrix ${\cal{L}}$, defined by
${\cal{L}}a_{\beta}=\lambda_{\beta}a_{\beta}$. Formally we can write the
solution of the master equation as
$|\rho(t)\rangle=\sum_{\beta}c_{\beta}\exp\left(-\lambda_{\beta}t\right)a_{\beta}$
where the coefficients $c_{\beta}$ are determined by the initial condition.
The eigenvalues $\lambda_{\beta}$ are real numbers since ${\cal{L}}$ is
symmetric and positive to have positive probabilities. Besides, one
eigenvalue, say $\beta=0$, must be equal to zero, $\lambda_{0}=0$, which
corresponds to the stationary case. The stationary probabilities can be
obtained from the left $\langle 0_{L}|$ and right $|0_{R}\rangle$ null
eigenvectors of ${\cal{L}}$ corresponding to zero eigenvalue
$|\rho^{(st)}\rangle=\frac{1}{\langle
0_{L}|0_{R}\rangle}|0_{R}\rangle=|0_{R}\rangle$.

Inasmuch as we are interested in the charge and heat noise for both dots we
introduce the particle and energy transport on the same footing, we need to
introduce both the particle-number and energy resolved density matrix

$\rho_{\gamma\eta}(\\{\mathbb{N},\mathbb{E}\\};t)\equiv\rho_{\gamma\eta}(\\{N_{r},N_{h},N_{c},E_{r},E_{h},E_{c}\\};t)$
(63)

which gives the probability of having a particle and an energy $N_{\alpha}$,
$E_{\alpha}$ in each corresponding lead $\alpha$ by the time $t$, where the
subscripts $\gamma$ and $\eta$ denote the different dot system states. In our
study we focus on the dynamics of the population, thus
$\left.\rho_{\gamma\eta}\right|_{\gamma=\eta}\equiv\rho_{\eta}$. $N_{r}$
($E_{r}$) distinguishes the charge (energy) at the reference reservoir for the
top quantum dot from the charge $N_{h,c}$ ($E_{h,c}$) at the hot and cold
reservoirs for the bottom quantum dot. Its Fourier transform is defined as

$\tilde{\rho}_{\eta}(\\{\chi,\kappa\\};t)=\sum_{\\{\mathbb{N},\mathbb{E}\\}}\rho_{\eta}(\\{\mathbb{N},\mathbb{E}\\};t)e^{i\chi\cdot\mathbb{N}}e^{i\kappa\cdot\mathbb{E}}$
(64)

where $\mathbb{N}=(N_{r},N_{h},N_{c})$ and $\mathbb{E}=(E_{r},E_{h},E_{c})$.
The counting fields are the conjugate particle and energy variables: for the
particle $\chi=(\chi_{r},\chi_{h},\chi_{c})$, and for the energy
$\mathbb{\kappa}=(\kappa_{r}$, $\kappa_{h},\kappa_{c})$ counting fields. After
Fourier transforming $\rho_{\eta}(\\{\mathbb{N},\mathbb{E}\\})$, the master
equation reads

$\frac{d}{dt}\tilde{\rho}\left(\\{\chi,\kappa\\};t\right)=-{\cal{L}}\left(\\{\chi,\kappa\\}\right)\tilde{\rho}\left(\\{\chi,\kappa\\};t\right)$
(65)

The explicit form for ${\cal{L}}\left(\\{\chi,\kappa\\}\right)$ is written in
Sec.VII) Here, the particle-energy number resolved density matrix fulfils the
master equation:
$\frac{d}{dt}\tilde{\rho}(\\{\chi,\kappa\\};t)=-{\cal{L}}(\\{\chi,\kappa\\})\tilde{\rho}(\\{\chi,\kappa\\};t)$.
From these equations the Liouvillian operator can be written as
${\cal{L}}(\\{\chi,\kappa\\})={\cal{L}}(\\{\mathbb{0},\mathbb{0}\\})+\tilde{{\cal{L}}}(\\{\chi,\kappa\\})$.
For the cumulant calculation we first compute the projection and pseudoinverse
operators ${\cal{P}}=1-|0_{R}\rangle\langle
0_{L}|,\quad{\cal{Q}}=1-{\cal{P}},\quad{\cal{R}}={\cal{Q}}{\cal{L}}^{-1}{\cal{Q}}$
where the Drazine inverse [53] of ${\cal{L}}$ is denoted as ${\cal{L}}^{-1}$.
For the top dot the charge and energy current at the reference lead it can be
derived from the first order cumulants. For instance, if one is interested in
computing charge current via lead $\alpha\in\\{r,h,c\\}$,

$\displaystyle j^{q}_{\alpha}=eC^{q}_{\alpha}=-e\langle
0_{L}|\left.\frac{\partial\tilde{{\cal{L}}}(\\{\chi,\kappa\\})}{\partial(i\chi_{\alpha})}\right|_{\chi=\kappa=\mathbb{0}}|0_{R}\rangle,$
(66)

which simply computes the change of particle number in lead $\alpha$. If one
is interested in computing the energy current via lead $\beta\in\\{r,h,c\\}$,

$\displaystyle j^{e}_{\beta}=C^{e}_{\beta}=-\langle
0_{L}|\left.\frac{\partial\tilde{{\cal{L}}}(\\{\chi,\kappa\\})}{\partial(i\kappa_{\beta})}\right|_{\chi=\kappa=\mathbb{0}}|0_{R}\rangle,$
(67)

which computes the change of total energy in lead $\beta$. The heat current
via lead $\gamma\in\\{r,h,c\\}$ is

$\displaystyle j^{h}_{\gamma}$
$\displaystyle=j^{e}_{\gamma}-\mu_{\alpha}j^{q}_{\gamma},$ (68)
$\displaystyle=C^{e}_{\gamma}-e\mu_{\gamma}C^{q}_{\gamma},$ (69)

Next, let us compute the charge and energy noise associated with lead $\alpha$
and $\beta$, respectively. For the evaluation of multidimensional TUR, one
must include the cross current-current correlation among three leads,
$\alpha\neq\beta$. First, the correlation of charge current via lead $\alpha$
and $\beta$ is

$\displaystyle C^{qq}_{\alpha\beta}$ $\displaystyle=-\langle
0_{L}|\left[\partial_{i\chi_{\alpha}}\partial_{i\chi_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle$
$\displaystyle+\left[\langle
0_{L}|\left[\partial_{i\chi_{\alpha}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}{\cal{R}}\left[\partial_{i\chi_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle+\left(\chi_{\alpha}\Leftrightarrow\chi_{\beta}\right)\right]$

where
$\left[\cdot\right]_{\mathbb{0}}=\left[\cdot\right]_{\chi=\kappa=\mathbb{0}}$.
When $\alpha=\beta$, Eq.(VI) provides the charge noise in lead $\alpha$. The
correlation of energy current via lead $\alpha$ and $\beta$ is similarly
written:

$\displaystyle C^{ee}_{\alpha\beta}$ $\displaystyle=-\langle
0_{L}|\left[\partial_{i\kappa_{\alpha}}\partial_{i\kappa_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle$
$\displaystyle+\left[\langle
0_{L}|\left[\partial_{i\kappa_{\alpha}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}{\cal{R}}\left[\partial_{i\kappa_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle+\left(\kappa_{\alpha}\Leftrightarrow\kappa_{\beta}\right)\right],$

where the counting field is simply replaced by $\kappa_{\alpha,\beta}$. The
correlation between the charge current via lead $\alpha$ and the energy
current via lead $\beta$ is as follows:

$\displaystyle C^{qe}_{\alpha\beta}$ $\displaystyle=-\langle
0_{L}|\left[\partial_{i\chi_{\alpha}}\partial_{i\kappa_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle$
$\displaystyle+\left[\langle
0_{L}|\left[\partial_{i\chi_{\alpha}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}{\cal{R}}\left[\partial_{i\kappa_{\beta}}\tilde{{\cal{L}}}\right]_{\mathbb{0}}|0_{R}\rangle+\left(\chi_{\alpha}\Leftrightarrow\kappa_{\beta}\right)\right].$

The evaluation of heat noise can be obtained from the charge and energy cross
correlations.

$\displaystyle C^{hh}_{\alpha\beta}$
$\displaystyle=C^{ee}_{\alpha\beta}-e\mu_{\beta}C^{eq}_{\alpha\beta}-e\mu_{\alpha}C^{qe}_{\alpha\beta}+e^{2}\mu_{\alpha}\mu_{\beta}C^{qq}_{\alpha\beta},$
(73)

which follows from the definition of heat current. The expression above is
enough for our purpose of computing the MTUR from heat currents via three
leads. Note that other combinations include

$\displaystyle C^{hq}_{\alpha\beta}$
$\displaystyle=C^{eq}_{\alpha\beta}-e\mu_{\alpha}C^{qq}_{\alpha\beta},$ (74)
$\displaystyle C^{qh}_{\alpha\beta}$
$\displaystyle=C^{qe}_{\alpha\beta}-e\mu_{\beta}C^{qq}_{\alpha\beta}.$ (75)

##  VII Quantum master equation

In this section we provide the derivation of the master equation employed in
the main text from the Lindblad equation. In the first subsection, We consider
the top quantum dot coupled to one metallic lead (reference reservoir) and one
superconducting lead. The master equation is derived for
$\Gamma_{r}\ll\Gamma_{S}$ so that the top quantum dot remains in the
superconducting phase. The dimension of Hilbert space is four,
$\ket{n_{t\uparrow}}\in\\{0,1\\},n_{t\downarrow}\in\\{0,1\\}$.

In the second subsection, the bottom quantum dot is introduced,
$n_{b}\in\\{0,1\\}$, by which the dimension of the Hilbert space is doubled.
There is one space with $n_{b}=0$, and the other space with $n_{b}=1$. The two
spaces are coupled by hot and cold lead which exchange electron and hole with
the bottom quantum dot. We work in the regime where
$T_{r,h,c}\gg\Gamma_{r,h,c}$, which allows us to neglect the off-diagonal
elements of the density matrix. Note that the master equation (rate equation
approach) is employed in Ref. [45], and here we provide its justification.

###  VII.1 A single quantum dot

In this section we obtain the quantum master equation for the simple case of a
single quantum dot attached to two contacts, one is a metallic reservoir and
the other is superconduting (N-QD-S). For such purpose we employ the quantum
Lindblard formalism . Our starting point is the effective Hamiltonian by the
superconducting proximity effect [40] for the N-QD-S system valid when the
largest energy scale $\Delta\rightarrow\infty$ as discussed in the main text.
Such Hamiltonian is

$\hat{H}_{n_{b}=0}=\varepsilon_{t}\sum_{\sigma=\uparrow,\downarrow}d^{\dagger}_{\sigma}d_{\sigma}+U_{t}\hat{n}_{\uparrow}\hat{n}_{\downarrow}+\Gamma_{S}d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}+\Gamma_{S}d_{\downarrow}d_{\uparrow},$
(76)

where the subscript $n_{b}=0$ is used because the Hamiltonian corresponds to
the situation where the particle number of the bottom quantum dot is fixed to
$n_{b}=0$. Writing this in matrix form in the electron number basis,
$|n_{\downarrow},n_{\uparrow}\rangle=\\{|0,0\rangle,|1,1\rangle,|1,0\rangle,|0,1\rangle\\}$,

$\displaystyle\hat{H}_{n_{b}=0}=\left(\begin{array}[]{cccc}0&\Gamma_{S}&0&0\\\
\Gamma_{S}&2\varepsilon_{t}+U_{t}&0&0\\\ 0&0&\varepsilon_{t}&0\\\
0&0&0&\varepsilon_{t}\end{array}\right).$ (81)

The eigenvalues of the first 2$\times$2 block are

$\varepsilon_{\pm}=\varepsilon_{t}+\frac{U_{t}}{2}\pm\sqrt{\Gamma_{S}^{2}+\left(\varepsilon_{t}+\frac{U_{t}}{2}\right)^{2}}.$
(82)

Corresponding eigenvectors are

$\ket{\pm}=u_{\pm}\ket{0,0}+v_{\pm}\ket{1,1},$ (83)

where

$u_{\pm}=\varepsilon_{\pm}/\sqrt{\Gamma_{S}^{2}+\varepsilon_{\pm}^{2}},\quad
v_{\pm}=-\Gamma_{S}/\sqrt{\Gamma_{S}^{2}+\varepsilon_{\pm}^{2}},$ (84)

where note that we used different sign of $v_{\pm}$ from the main text to
simplify the notation in the following discussion. The master equation
Eq.(123) that we obtain later is independent of this choice. Introducing the
unitary matrix

$\displaystyle U_{0}=\left(\begin{array}[]{cccc}u_{+}&u_{-}&0&0\\\
v_{+}&v_{-}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{array}\right),\,\,U_{0}^{\dagger}=\left(\begin{array}[]{cccc}u^{*}_{+}&v^{*}_{+}&0&0\\\
u^{*}_{-}&v^{*}_{-}&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{array}\right),$ (93)

the Hamiltonian can be diagonalized:
$D=U_{0}^{\dagger}H_{n_{b}=0}U_{0}=\text{diag}(\varepsilon_{+},\varepsilon_{-},\varepsilon_{t},\varepsilon_{t})$.
Let us explicitly write the operators in the basis of
$\ket{n_{\downarrow},n_{\uparrow}}$:

$\displaystyle\hat{d}_{\uparrow}=\left(\begin{array}[]{cccc}0&0&0&1\\\
0&0&0&0\\\ 0&1&0&0\\\
0&0&0&0\end{array}\right),\,\,\,\,\hat{d}^{\dagger}_{\uparrow}=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&0&1&0\\\ 0&0&0&0\\\ 1&0&0&0\end{array}\right),$ (102)
$\displaystyle\hat{d}_{\downarrow}=\left(\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&0\\\ 0&0&0&0\\\
0&1&0&0\end{array}\right),\,\,\,\,\hat{d}^{\dagger}_{\downarrow}=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&0&0&1\\\ 1&0&0&0\\\ 0&0&0&0\end{array}\right).$ (111)

which yields straightforward expressions for number operators:
$d^{\dagger}_{\uparrow}d_{\uparrow}=\text{diag}(0,1,0,1)$, and
$d^{\dagger}_{\downarrow}d_{\downarrow}=\text{diag}(0,1,1,0)$. The coupling to
the normal lead induces the transition between the eigenstates. Since the
normal lead is eligible to exchange an electron, to construct the Lindblad
operator we need to express $d_{\sigma}^{(\dagger)}$ in terms of eigen
operators of the system. That is,
$\tilde{d}_{\uparrow}=U_{0}^{\dagger}d_{\uparrow}U_{0}$,
$\tilde{d}_{\downarrow}=U_{0}^{\dagger}d_{\downarrow}U_{0}$. The Lindblad
operator comes with the coupling strength between the system and the
electronic bath,

$\sqrt{\gamma^{e}_{N,i\leftarrow j}}=\sqrt{\Gamma_{r}f_{r}^{e}(E_{i}-E_{j})},$
(112)

where $f^{e}_{r}(E_{i}-E_{j})=[e^{(E_{i}-E_{j}-\mu_{r})/T_{r}}+1]^{-1}$ is the
Fermi-Dirac distribution, for a transition of adding an electron to the system
that accompanies energy change $(E_{i}-E_{j})$ in the system. For example, the
Lindblad operator that removes an spin-up electron from the system to the
reservoir is:

$\displaystyle\tilde{L}_{\uparrow}=\left(\begin{array}[]{cccc}0&0&0&u_{+}^{*}\sqrt{\gamma^{h}_{+\uparrow}}\\\
0&0&0&u_{-}^{*}\sqrt{\gamma^{h}_{-\downarrow}}\\\
v_{+}\sqrt{\gamma^{h}_{\downarrow+}}&v_{-}\sqrt{\gamma^{h}_{\downarrow-}}&0&0\\\
0&0&0&0\end{array}\right)$ (117)

where
$\gamma^{e(h)}_{ij}\equiv\gamma_{r,|i\rangle\leftarrow|j\rangle}=\Gamma_{r}f_{r}^{e(h)}(E_{i}-E_{j})$
($e$ for electrons and $h$ for holes) and
$i,j\in\\{+,-,\downarrow,\uparrow\\}$. Other Lindblad operators
$\tilde{L}^{\dagger}_{\uparrow},\tilde{L}_{\downarrow},\tilde{L}^{\dagger}_{\downarrow}$
are expressed similarly. The equation of motion of the density matrix is

$\frac{d\rho}{dt}=-i[H_{\mathrm{top}},\rho]+\mathcal{D}[L_{\sigma}]+\mathcal{D}[L^{\dagger}_{\sigma}],$
(118)

where the dissipator operator due to the normal lead is

$\mathcal{D}[L^{\dagger}_{\sigma}]=L_{\sigma}\rho
L^{\dagger}_{\sigma}-\frac{1}{2}\\{L^{\dagger}_{\sigma}L_{\sigma},\rho\\}.$
(119)

Rewriting the equation of motion in the eigen state basis of Hamiltonian,

$\displaystyle U_{0}^{\dagger}\frac{d\rho}{dt}U_{0}$
$\displaystyle=U_{0}^{\dagger}\left(-i[H_{\mathrm{top}},\rho]+\mathcal{D}[L_{\sigma}]+\mathcal{D}[L^{\dagger}_{\sigma}]\right)U_{0},$
(120) $\displaystyle\frac{d\tilde{\rho}}{dt}$
$\displaystyle=-i[D,\tilde{\rho}]+\tilde{\mathcal{D}}[L_{\sigma}]+\tilde{\mathcal{D}}[L^{\dagger}_{\sigma}],$
(121)

where

$\displaystyle\tilde{\mathcal{D}}[L^{\dagger}_{\sigma}]$
$\displaystyle=U_{0}^{\dagger}{\mathcal{D}}[L^{\dagger}_{\sigma}]U_{0}=(U_{0}^{\dagger}L_{\sigma}U_{0})\tilde{\rho}(U_{0}^{\dagger}L^{\dagger}_{\sigma}U_{0})$
$\displaystyle-\frac{1}{2}\\{U_{0}^{\dagger}(L^{\dagger}_{\sigma}L_{\sigma})U_{0},\tilde{\rho}\\}.$
(122)

For the derivation of the master equation from the quantum Lindblad equation
[Eq. (118)] focusing on the change of the occupation numbers of eigen states
and disregarding the coherences. For that issue we employ the density matrix
with diagonal elements only $\tilde{\rho}=U_{0}^{\dagger}\rho
U_{0}=\text{diag}(\rho_{+},\rho_{-},\rho_{\downarrow},\rho_{\uparrow})$. This
approach is justified as long as $\Gamma_{S}\ll k_{B}T_{\alpha}$ ($\alpha\in
r,h,c$). Considering this approach we obtain the following set of equation of
motions:

$\displaystyle\frac{\partial\rho_{+}}{\partial
t}=\left[|u_{+}|^{2}(\gamma^{h}_{+,\uparrow}\rho_{\uparrow}+\gamma^{h}_{+,\downarrow}\rho_{\downarrow})-|v_{+}|^{2}(\gamma^{h}_{\downarrow+}+\gamma^{h}_{\uparrow+})\rho_{+}\right]+\left[|v_{+}|^{2}(\gamma^{e}_{+,\uparrow}\rho_{\uparrow}+\gamma^{e}_{+,\downarrow}\rho_{\downarrow})-|u_{+}|^{2}(\gamma^{e}_{\downarrow+}+f^{e}_{\uparrow+})\rho_{+}\right],$
(123) $\displaystyle\frac{\partial\rho_{-}}{\partial
t}=\left[|u_{-}|^{2}(\gamma^{h}_{-\uparrow}\rho_{\uparrow}+\gamma^{h}_{-\downarrow}\rho_{\downarrow})-|v_{-}|^{2}(\gamma^{h}_{\downarrow-}+\gamma^{h}_{\uparrow-})\rho_{-}\right]+\left[|v_{-}|^{2}(\gamma^{e}_{-\uparrow}\rho_{\uparrow}+\gamma^{e}_{-\downarrow}\rho_{\downarrow})-|u_{-}|^{2}(\gamma^{e}_{\downarrow-}+\gamma^{e}_{\uparrow-})\rho_{-}\right],$
(124) $\displaystyle\frac{\partial\rho_{\downarrow}}{\partial
t}=\left[|v_{+}|^{2}\gamma^{h}_{\downarrow+}\rho_{+}+|v_{-}|^{2}\gamma^{h}_{\downarrow-}\rho_{-}-(|u_{+}|^{2}\gamma^{h}_{+\downarrow}+|u_{-}|^{2}\gamma^{h}_{-\downarrow})\rho_{\downarrow}\right]+[|u_{+}|^{2}\gamma^{e}_{\downarrow+}\rho_{+}+|u_{-}|^{2}\gamma^{e}_{\downarrow-}\rho_{-}-(|v_{+}|^{2}\gamma^{e}_{+\downarrow}+|v_{-}|^{2}\gamma^{e}_{-\downarrow})\rho_{\downarrow}],$
(125) $\displaystyle\frac{\partial\rho_{\uparrow}}{\partial
t}=\left[|v_{+}|^{2}\gamma^{h}_{\uparrow+}\rho_{+}+|v_{-}|^{2}\gamma^{h}_{\uparrow-}\rho_{-}-(|u_{+}|^{2}\gamma^{h}_{+\uparrow}+|u_{-}|^{2}\gamma^{h}_{-\uparrow})\rho_{\uparrow}\right]+[|u_{+}|^{2}\gamma^{e}_{\uparrow+}\rho_{+}+|u_{-}|^{2}\gamma^{e}_{\uparrow-}\rho_{-}-(|v_{+}|^{2}\gamma^{e}_{+\uparrow}+|v_{-}|^{2}\gamma^{e}_{-\uparrow})\rho_{\uparrow}].$
(126)

where the amplitude of each term can be read from
$\ket{\pm}=u_{\pm}\ket{0,0}+v_{\pm}\ket{1,1}$. When single electron states
$\ket{\uparrow}$ or $\ket{\downarrow}$ receives one electron to become
$\ket{\pm}$, the transition rate is $|v_{\pm}|^{2}$. When single electron
states $\ket{\uparrow}$ or $\ket{\downarrow}$ emit one electron to become
$\ket{\pm}$, the transition rate is $|u_{\pm}|^{2}$. This explains the first
two equations with positive signs. The other terms can be read similarly.

###  VII.2 A double quantum dot

The extension of the previous derivation of the master equation for the double
quantum dot system is straightforward. We simply need to extend the Hilbert
space to take into account the occupancy of the bottom quantum dot,
$n_{b}\in\\{0,1\\}$. Our Hamiltonian is
$H_{DQD}=\text{blkdiag}(H_{n_{b}=0},H_{n_{b}=1})$. The first block diagonal
Hamiltonian $H_{n_{b}=0}$ is shown in Eq. (81). The second block diagonal
Hamiltonian is the following:

$\displaystyle
H_{\mathrm{n_{b}=1}}=\left(\begin{array}[]{cccc}0&\Gamma_{S}&0&0\\\
\Gamma_{S}&2\varepsilon_{t}+U_{t}+2U_{tb}&0&0\\\
0&0&\varepsilon_{t}+\varepsilon_{b}+U_{tb}&0\\\
0&0&0&\varepsilon_{t}+\varepsilon_{b}+U_{tb}\end{array}\right).$ (131)

with eigenvalues
$\varepsilon_{\pm}^{\prime}=\varepsilon_{t}+U_{tb}+\frac{U_{t}}{2}\pm\sqrt{\Gamma_{S}^{2}+(\varepsilon_{t}+U_{tb}+\frac{U_{t}}{2})^{2}}$.
The corresponding eigenvectors are
$\ket{\pm}=u_{\pm}\ket{0,0}+v_{\pm}\ket{1,1}$. By introducing the unitary
matrix,

$\displaystyle U_{1}=\left(\begin{array}[]{cccc}u_{+,1}&u_{-,1}&0&0\\\
v_{+,1}&v_{-,1}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{array}\right),\,\,\,\,U_{1}^{\dagger}=\left(\begin{array}[]{cccc}u^{*}_{+,1}&v^{*}_{+,1}&0&0\\\
u^{*}_{-,1}&v^{*}_{-,1}&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{array}\right),$ (140)

we have the diagonalization of the Hamiltonian:
$D_{n_{b}=1}=U_{1}^{\dagger}H_{n_{b}=1}U_{1}=\text{diag}(\varepsilon_{+},\varepsilon_{-},\varepsilon_{t}+\varepsilon_{b}+U_{tb},\varepsilon_{t}+\varepsilon_{b}+U_{tb})$.
The Hamiltonian of the double quantum dot is digonalized by the unitary matrix
$U=\text{blkdiag}(U_{0},U_{1})$. In the eigen state basis, the Hamiltonian is
diagonal: $D_{tb}=U^{\dagger}H_{tb}U=\text{blkdiag}(D_{n_{b}=0},D_{n_{b}=1})$.
We then need to write the equation of motion for the density matrix ignoring
the coherence terms
$\rho=\text{diag}(\rho_{+,0},\rho_{-,0},\rho_{\downarrow,0},\rho_{\uparrow,0},\rho_{+,1},\rho_{-,1},\rho_{\downarrow,1},\rho_{\uparrow,1})$.
Because the 8$\times$8 unitary matrix is block diagonal, our previous
expressions for the single quantum dot will equally apply for sector $n_{b}=0$
and $n_{b}=1$ with modified eigenvalues and eigenvectors. Terms newly added
are the creation and removal of electrons in the bottom quantum dot. The
creation and annihilation operators of an electron in the bottom quantum dot
are expressed as follows

$\displaystyle d^{\dagger}_{b}$ $\displaystyle=\left(\begin{array}[]{c}0\\\
\mathcal{I}_{4}\end{array}\right),\,\,\,\,\,\,d_{b}=\left(\begin{array}[]{c}\mathcal{I}_{4}\\\
0\end{array}\right),$ (145) $\displaystyle\tilde{d}^{\dagger}_{b}$
$\displaystyle=U^{\dagger}d^{\dagger}_{t}U=\left(\begin{array}[]{cc}0&0\\\
U_{1}^{\dagger}U_{0}&0\end{array}\right),$ (148) $\displaystyle\tilde{d}_{b}$
$\displaystyle=U^{\dagger}d_{b}U=\left(\begin{array}[]{cc}0&U_{0}^{\dagger}U_{1}\\\
0&0\end{array}\right),$ (151)

where the first two are in the electron number basis
$\ket{n_{\downarrow},n_{\uparrow},n_{b}}$, and the latter two are the
expressions in the eigen basis. From these, we can build the explicit
expression of dissipators from the left and right lead coupled to the bottom
quantum dot

$\displaystyle\tilde{\mathcal{D}}[d^{\dagger}_{b}]$
$\displaystyle=U^{\dagger}{\mathcal{D}}[d^{\dagger}_{b}]U=\tilde{d}_{b}\tilde{\rho}\tilde{d}^{\dagger}_{b}-\frac{1}{2}\\{\tilde{d}_{b}^{\dagger}\tilde{d}_{b},\tilde{\rho}\\}=$
(154)
$\displaystyle\left(\begin{array}[]{cc}U_{0}^{\dagger}U_{1}\rho_{n_{b}=1}U_{1}^{\dagger}U_{0}&0\\\
0&-\rho_{n_{b}=1}\end{array}\right),$

which corresponds to the physical process of removing one electron from the
bottom quantum dot. Thus, it is related to the density matrix sector of
$(n_{b}=1)$. The other dissipator is to add a particle to the bottom quantum
dot

$\displaystyle\tilde{\mathcal{D}}[d_{b}]$
$\displaystyle=U^{\dagger}{\mathcal{D}}[d_{b}]U=\tilde{d}^{\dagger}_{b}\tilde{\rho}\tilde{d}_{b}-\frac{1}{2}\\{\tilde{d}_{b}\tilde{d}^{\dagger}_{b},\tilde{\rho}\\}=$
(157) $\displaystyle\left(\begin{array}[]{cc}-\rho_{n_{b}=0}&0\\\
0&U_{1}^{\dagger}U_{0}\rho_{n_{b}=0}U_{0}^{\dagger}U_{1}\end{array}\right),$

From the two dissipators, there are scattering between states $\ket{+,1}$ and
$\ket{-,0}$, for example. Specifically, the diagonal elements are the
following:

$\displaystyle\left(U_{0}^{\dagger}U_{1}\rho_{n_{b}=1}U_{1}^{\dagger}U_{0}\right)_{11}$
$\displaystyle=|u_{+,0}^{*}u_{+,1}+v_{+,0}^{*}v_{+,1}|^{2}\rho_{+,1}+|u_{+,0}^{*}u_{-,1}+v_{+,0}^{*}v_{-,1}|^{2}\rho_{-,1},$
(158)
$\displaystyle\left(U_{0}^{\dagger}U_{1}\rho_{n_{b}=1}U_{1}^{\dagger}U_{0}\right)_{22}$
$\displaystyle=|u_{-,0}^{*}u_{-,1}+v_{-,0}^{*}v_{-,1}|^{2}\rho_{-,1}+|u_{-,0}^{*}u_{+,1}+v_{-,0}^{*}v_{+,1}|^{2}\rho_{+,1},$
(159)
$\displaystyle\left(U_{1}^{\dagger}U_{0}\rho_{n_{b}=0}U_{0}^{\dagger}U_{1}\right)_{11}$
$\displaystyle=|u_{+,1}^{*}u_{+,0}+v_{+,1}^{*}v_{+,0}|^{2}\rho_{+,0}+|u_{+,1}^{*}u_{-,0}+v_{+,1}^{*}v_{-,0}|^{2}\rho_{-,0},$
(160)
$\displaystyle\left(U_{1}^{\dagger}U_{0}\rho_{n_{b}=0}U_{0}^{\dagger}U_{1}\right)_{22}$
$\displaystyle=|u_{-,1}^{*}u_{-,0}+v_{-,1}^{*}v_{-,0}|^{2}\rho_{-,0}+|u_{-,1}^{*}u_{+,0}+v_{-,1}^{*}v_{+,0}|^{2}\rho_{+,0},$
(161)

For different scatterings, the associated energy transfer is different. That
is, we need to introduce different argument for the Fermi-Dirac distribution
that is multiplied for each scattering process. Thus, we need to include the
coupling strength between the system and baths
$\sim\sqrt{\Gamma_{l}f^{e,h}_{l,\ket{i}\leftarrow\ket{j}}}$ to complete the
expression of the dissipators. As a result, for the double quantum dot system
we have additional terms accounting for the scatterings changing the electron
number of bottom quantum dot, $n_{b}$.

$\displaystyle\frac{\partial\rho_{+,0}}{\partial
t}=\left.\frac{\partial\rho_{+,0}}{\partial t}\right|_{n_{b}=0\leftarrow 0}$
$\displaystyle+\sum_{l=h,c}\gamma^{h}_{l,\ket{+,0}\leftarrow\ket{+,1}}|u_{+,0}^{*}u_{+,1}+v_{+,0}^{*}v_{+,1}|^{2}\rho_{+,1}+\gamma^{h}_{l,\ket{+,0}\leftarrow\ket{-,1}}|u_{+,0}^{*}u_{-,1}+v_{+,0}^{*}v_{-,1}|^{2}\rho_{-,1}$
$\displaystyle-\sum_{l=h,c}\gamma^{e}_{l,\ket{+,1}\leftarrow\ket{+,0}}|u_{+,1}^{*}u_{+,0}+v_{+,1}^{*}v_{+,0}|^{2}\rho_{+,0}+\gamma^{e}_{l,\ket{-,1}\leftarrow\ket{+,0}}|u_{-,1}^{*}u_{+,0}+v_{-,1}^{*}v_{+,0}|^{2}\rho_{+,0},$
$\displaystyle\frac{\partial\rho_{-,0}}{\partial
t}=\left.\frac{\partial\rho_{-,0}}{\partial t}\right|_{n_{b}=0\leftarrow 0}$
$\displaystyle+\sum_{l=h,c}\gamma^{h}_{l,\ket{-,0}\leftarrow\ket{-,1}}|u_{-,0}^{*}u_{-,1}+v_{-,0}^{*}v_{-,1}|^{2}\rho_{-,1}+\gamma^{h}_{l,\ket{-,0}\leftarrow\ket{+,1}})|u_{-,0}^{*}u_{+,1}+v_{-,0}^{*}v_{+,1}|^{2}\rho_{+,1}$
$\displaystyle-\sum_{l=h,c}\gamma^{e}_{l,\ket{-,1}\leftarrow\ket{-,0}}|u_{-,1}^{*}u_{-,0}+v_{-,1}^{*}v_{-,0}|^{2}\rho_{-,0}+\gamma^{e}_{l,\ket{+,1}\leftarrow\ket{-,0}}|u_{+,1}^{*}u_{-,0}+v_{+,1}^{*}v_{-,0}|^{2}\rho_{-,0},$
$\displaystyle\frac{\partial\rho_{\downarrow,0}}{\partial
t}=\left.\frac{\partial\rho_{\downarrow,0}}{\partial
t}\right|_{n_{b}=0\leftarrow 0}$
$\displaystyle+\sum_{l=h,c}\gamma^{h}_{l,\ket{\downarrow,0}\leftarrow\ket{\downarrow,1}}\rho_{\downarrow,1}-\gamma^{e}_{l,\ket{\downarrow,1}\leftarrow\ket{\downarrow,0}}\rho_{\downarrow,0},$
$\displaystyle\frac{\partial\rho_{\uparrow,0}}{\partial
t}=\left.\frac{\partial\rho_{\uparrow,0}}{\partial
t}\right|_{n_{b}=0\leftarrow 0}$
$\displaystyle+\sum_{l=h,c}\gamma^{h}_{l,\ket{\uparrow,0}\leftarrow\ket{\uparrow,1}}\rho_{\uparrow,1}-\gamma^{e}_{l,\ket{\uparrow,1}\leftarrow\ket{\uparrow,0}}\rho_{\uparrow,0},$

where $\left.\frac{\partial\rho_{k}}{\partial t}\right|_{n_{b}=0\leftarrow 0}$
are the same equations as appears in Eq. (123-126) with $n_{b}=0$ in which the
transition rates read as follows:
$\gamma^{e(h)}_{ij}=\gamma^{e(h)}_{|i,n_{b}\rangle\leftarrow|jn_{b}\rangle}=\Gamma_{r}f^{e(h)}(E_{i,n_{b}}-E_{j,n_{b}})$
with $i,j\in\\{\pm,\sigma\\}$, the particle number $n_{b}$ in the bottom dot
does not change. We have another set of equation of motion for $n_{b}=1$:

$\displaystyle\frac{\partial\rho_{+,1}}{\partial
t}=\left.\frac{\partial\rho_{+,1}}{\partial t}\right|_{n_{b}=1\leftarrow 1}$
$\displaystyle+\sum_{l=h,c}\gamma^{e}_{l,\ket{+,1}\leftarrow\ket{+,0}}|u_{+}^{{}^{\prime}*}u_{+}+v_{+}^{{}^{\prime}*}v_{+}|^{2}\rho_{+}+\gamma^{e}_{l,\ket{+,1}\leftarrow\ket{-,0}}|u_{+,1}^{*}u_{-,0}+v_{+,1}^{*}v_{-,0}|^{2}\rho_{-,0}$
$\displaystyle-\sum_{l=h,c}\gamma^{h}_{l,\ket{+,0}\leftarrow\ket{+,1}}|u_{+,0}^{*}u_{+,1}+v_{+,0}^{*}v_{+,1}|^{2}\rho_{+,1}+\gamma^{h}_{l,\ket{-,0}\leftarrow\ket{+,1}}|u_{-,0}^{*}u_{+,1}+v_{-,0}^{*}v_{,1}+|^{2}\rho_{+,1},$
$\displaystyle\frac{\partial\rho_{-,1}}{\partial
t}=\left.\frac{\partial\rho_{-,1}}{\partial t}\right|_{n_{b}=1\leftarrow 1}$
$\displaystyle+\sum_{l=h,c}\gamma^{e}_{l,\ket{-,1}\leftarrow\ket{-,0}}|u_{-,1}^{*}u_{-,0}+v_{-,1}^{*}v_{-,0}|^{2}\rho_{-,0}+\gamma^{e}_{l,\ket{-,1}\leftarrow\ket{+,0}}|u_{-,1}^{*}u_{+,0}+v_{-,1}^{*}v_{+,0}|^{2}\rho_{+,0}$
$\displaystyle-\sum_{l=h,c}\gamma^{h}_{l,\ket{-,0}\leftarrow\ket{-,1}}|u_{-,0}^{*}u_{-,1}+v_{-,0}^{*}v_{-,1}|^{2}\rho_{-,1}+\gamma^{h}_{l,\ket{+,0}\leftarrow\ket{-,1}}|u_{+,0}^{*}u_{-,1}+v_{+,0}^{*}v_{-,1}|^{2}\rho_{-,1},$
$\displaystyle\frac{\partial\rho_{\downarrow,0}}{\partial
t}=\left.\frac{\partial\rho_{\downarrow,1}}{\partial
t}\right|_{n_{b}=1\leftarrow 1}$
$\displaystyle+\sum_{l=h,c}\gamma^{e}_{l,\ket{\downarrow,1}\leftarrow\ket{\downarrow,0}}\rho_{\downarrow,0}-\gamma_{l,\ket{\downarrow,0}\leftarrow\ket{\downarrow,1}}^{h}\rho_{\downarrow,1},$
$\displaystyle\frac{\partial\rho_{\uparrow,1}}{\partial
t}=\left.\frac{\partial\rho_{\uparrow,1}}{\partial
t}\right|_{n_{b}=1\leftarrow 1}$
$\displaystyle+\sum_{l=h,c}\gamma^{e}_{l,\ket{\uparrow,1}\leftarrow\ket{\uparrow,0}}\rho_{\uparrow,0}-\gamma^{h}_{l,\ket{\uparrow,0}\leftarrow\ket{\uparrow,1}}\rho_{\uparrow,1},$

where $\left.\frac{\partial\rho_{k}}{\partial t}\right|_{n_{b}=1\leftarrow 1}$
are the same equations as appears in Eq. (123-126) with $n_{b}=1$. As a
result, we obtain the master equation used in the main text.