††thanks<EMAIL_ADDRESS>
# Graph Theoretic Analysis of Three-Terminal Quantum Dot Thermocouples:
Onsager Relations and Spin-Thermoelectric Effects

Nikhil Gupt  Indian Institute of Technology Kanpur, Kanpur, Uttar Pradesh
208016, India  Shuvadip Ghosh  Indian Institute of Technology Kanpur, Kanpur,
Uttar Pradesh 208016, India  Arnab Ghosh  Indian Institute of Technology
Kanpur, Kanpur, Uttar Pradesh 208016, India

###### Abstract

We introduce a simplified model for a three-terminal quantum thermocouple
consisting of two strongly-coupled quantum dots. To elucidate spin-dependent
Seebeck and Peltier effects, we employ a microscopic Hamiltonian and map the
Lindblad master equation onto a quantum transition network, capturing the key
working principles for both reciprocal effects. Our analysis reveals quantum
thermodynamic networks encompassing both Coulomb interaction and spin-flipping
processes, lead to the emergence of spin-thermolectric effects. Using
algebraic graph theory, we recover the phenomenological law of irreversible
thermodynamics from the stochastic version of the entropy production rate
expressed in terms of cycle flux and cycle forces. Remarkably, Onsager
reciprocity and Kelvin relation for transport coefficients find their premises
in the properties of cycle flux trajectories within the quantum transition
network. This underscores the universal generality of thermodynamic principles
across classical and quantum realms, despite their fundamentally different
basis from classical laws of irreversible thermodynamics relying on local
equilibrium assumptions.

##  I Introduction

Thermoelectric devices have garnered significant attention owing to the
continual demand for innovative and effective approaches to temperature
sensors, heat pumps, and energy conversion Rowe (1995); DiSalvo (1999);
Goldsmid (2009); Shakouri and Zebarjadi (2009); Dubi and Di Ventra (2011);
Mazza  _et al._ (2015); Benenti  _et al._ (2017). This interest is rooted in
the phenomenon of thermoelectricity, where a temperature gradient induces an
electric current (Seebeck effect), and a potential gradient induces a heat
current (Peltier effect). From a thermodynamic point of view, a non-
equilibrium system experiences a distinct set of generalized thermodynamic
forces, arising from its simultaneous couplings with different reservoirs
Callen (1985); Landi and Paternostro (2021). The system’s response to these
external thermodynamic forces is reflected in a corresponding set of
generalized thermodynamic fluxes. The concept has been well investigated in
classical irreversible thermodynamics, with Onsager’s groundbreaking work on
the reciprocity principle of thermoelectric phenomena Onsager (1931a, b);
Callen (1948). Traditionally, thermocouples consisting of two different metal
wires, are used to observe such reciprocal effects. Only in recent times,
experimental research on magnetic metals and insulators, have experienced the
emergence of the spin Seebeck effect (SSE), wherein a spin current is
generated in response to a thermal gradient Uchida  _et al._ (2008); Wu  _et
al._ (2015); Zhou  _et al._ (2021), and conversely, the spin Peltier effect
(SPE), involves a spin voltage producing a thermal current Flipse  _et al._
(2014); Daimon  _et al._ (2016); Ohnuma  _et al._ (2017). The above findings
have ignited renewed enthusiasm among researchers to grasp the fundamental
aspects of spin caloritronics Bauer  _et al._ (2012); Boona  _et al._ (2014);
Ronetti  _et al._ (2016); ichi UCHIDA (2021) and explore practical
applications such as waste heat recovery and on-chip refrigeration for future
nanoelectronics. A a result, there is a considerable interest in understanding
the quantum thermodynamics of nanoscale thermoelectrics through theoretical
modelings Di Ventra (2008); Nazarov and Blanter (2009); Ihn (2010); Heikkilä
(2013); Ren (2013); Whitney  _et al._ (2016, 2018); Wang  _et al._ (2022) and
experimental setups involving quantum dot (QD) nanostructures, nanowires, and
two-dimensional materials van Houten  _et al._ (1992); Lee  _et al._ (2016);
Svilans  _et al._ (2016); Erlingsson  _et al._ (2017); Patel  _et al._ (2020);
Han  _et al._ (2020); Yang  _et al._ (2023).

The quantized energy levels and strong on-site Coulomb interactions among QDs,
make them excellent candidates for thermoelectric applications Esposito  _et
al._ (2009); Nakpathomkun  _et al._ (2010); Donsa  _et al._ (2014); Sothmann
_et al._ (2014); Whitney  _et al._ (2016); Erdman  _et al._ (2017); Whitney
_et al._ (2018); Wang  _et al._ (2022) and various other nanoscale thermal
devices Esposito  _et al._ (2012); Thierschmann  _et al._ (2015); Jiang  _et
al._ (2015); Zhang  _et al._ (2017); Ghosh  _et al._ (2022). While the
discrete QD spectrum can be fine-tuned via external gate voltages and offers
energy-selective transport, the strong Coulombic interaction between electrons
on capacitively coupled QDs can facilitate the transfer of precise amounts of
energy from the heat reservoirs. However, the use of QDs as working substances
for quantum thermodynamic devices, characterized by a limited number of
quantum states, necessitates a completely new understanding of these devices
Whitney  _et al._ (2018). The typical thermalization length being larger than
the nanoscale dimension forces these systems to behave in a highly non-trivial
manner, and their transport properties cannot be adequately described by the
usual Boltzmann transport equation Datta (2005), which primarily relies on the
local equilibrium assumptions.

On the contrary, the Lindblad master equation, formulated in terms of the
density matrix, is used as the preferred tool for examining the thermodynamic
properties of the open quantum systems Breuer and Petruccione (2007);
Gelbwaser-Klimovsky  _et al._ (2015); Joulain  _et al._ (2016); Ghosh  _et
al._ (2017); Potts (2019); Gupt  _et al._ (2021, 2022); Ghosh  _et al._
(2022). Though it is quite effective in accurately calculating the steady-
state currents amid non-equilibrium conditions, it does not reveal any
information about the operational principles and the nature of the transport
coefficients involved in complex quantum systems. In contrast, network theory
in recent years has emerged as a powerful instrument for comprehending non-
equilibrium quantum systems Wang  _et al._ (2022); Gupt  _et al._ (2023). In
this framework, dissipative quantum dynamics can be represented as a weighted
network featuring nodes and edges Schnakenberg (1976). Here, vertices (nodes)
signify quantum states, and edges denote non-equilibrium transitions from one
quantum state to another, with positive flux rates. Network theory has been
applied for many years to explore complex biological phenomena and chemical
reactions Hill and Chen (1975); Kohler and Vollmerhaus (1980); Ren (2017);
Dutta  _et al._ (2020). However, recent work by Wang et al. Wang  _et al._
(2022) has drawn huge attention by utilizing network theory to understand the
principle working mechanism of quantum thermal devices. The present authors
have extended the technique further to molecular systems to unravel hidden
electron transfer pathways in solar cells under strong non-equilibrium
conditions Gupt  _et al._ (2023).

In this paper, we leverage the advantages of network theory to elucidate the
operational principles of spin-thermoelectric effects within a three-terminal
quantum setup, closely resembling classical thermocouples. We demonstrate how
spin and energy currents, obtained from the quantum master equation, are
linked to the thermodynamic forces, manifesting spin-Seebeck and spin-Peltier,
as thermodynamic cross-effects. Close parallelisms between the microscopic and
the macroscopic description of the non-equilibrium system are established via
cycle force and cycle fluxes within a basic graph and thermodynamic forces and
fluxes of phenomenological laws. The central concept being used here is an
expression of the entropy production rate within the framework of the
algebraic graph.

The present work is organized as follows: In Sec. II, we introduce the basic
model of the quantum thermocouple and present the microscopic description
using the Lindblad master equation and quantum kinetic Pauli master equation.
We elaborate the basic framework of network theory in the context of spin-
thermoelectric effects in Sec. III and recover the phenomenological law of
irreversible thermodynamics and Onsager’s reciprocity in terms of network
cycle flux and forces. Operational principles of both spin-Seebeck and spin-
Peltier effects are presented in Sec. IV and finally, we conclude in Sec. V.

##  II Microscopic Model and Quantum Master equation

The basic model of a quantum thermocouple consists of two strongly coupled
quantum dots (QDs) via Coulomb interaction. The lower quantum dot, denoted as
${\rm QD}_{l}$, is simultaneously coupled with a spinful free-electron
reservoir (on the left) and a magnon bath (on the right), both maintained at
an equal temperature ($T_{0}$), as depicted in Fig. 1. The upper quantum dot,
${\rm QD}_{u}$, is only coupled with a spinless free-electron reservoir,
acting as a junction like in a classical thermocouple. The spinful free-
electron reservoir comprises spin-polarized electrons with both spin-up
($\uparrow$) and spin-down ($\downarrow$) orientations Vandaele  _et al._
(2017); Wang  _et al._ (2022). In contrast, the spinless free-electron
reservoir in the middle consists of electrons without any distinct spin
orientation, and the magnon bath at the right is responsible for inducing
spin-flipping of the ${\rm QD}_{l}$ electrons Ren (2013); Wang  _et al._
(2022); Vandaele  _et al._ (2017); Sothmann and Büttiker (2012). The three-
terminal QD model presented here bears a striking similarity to a classical
thermocouple, particularly in the manifestation of both the spin-Seebeck and
spin-Peltier effects (SSE and SPE). In SSE, a spin current emerges under the
influence of a temperature gradient ($\delta T$), while SPE occurs with the
application of a spin bias voltage at the lower terminals, resembling the open
ends of a conventional thermocouple. Notably, the difference between in the
statistical properties of the magnon and electron reservoirs plays a crucial
role in generating spin-thermoelectric effects. This distinction can be
likened to the role of dissimilar metal wires in a classical thermocouple,
highlighting its significance within this quantum framework.

Figure 1: Inset: Schematic diagram of the classical thermocouple. Main:
Schematic diagram of a three-terminal Coulomb-coupled QD thermocouple. The
lower quantum dot (${\rm QD}_{l}$) is coupled to the left reservoir i.e. a
spinful reservoir (in green) and the right reservoir (in blue) i.e. a magnon
bath. Both terminals are kept at equal temperatures and serve as cold ends.
The upper quantum dot (${\rm QD}_{u}$) is coupled to the middle reservoir (in
red) i.e. a spinless electron reservoir. Here, heat is transferred from the
middle reservoir which acts as a junction (hot end) and the spin current is
across the lower two terminals, analogous to the open ends of a thermocouple.

The total Hamiltonian of the entire three-terminal setup is given below

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{\rm S}+H_{\rm B}+H_{\rm
I},$ $\displaystyle H_{\rm S}$ $\displaystyle=$
$\displaystyle\sum_{\sigma=\\{\uparrow,\downarrow\\}}\varepsilon_{l\sigma}{n}_{l\sigma}+\varepsilon_{u}{n}_{u}+\sum_{\sigma=\\{\uparrow,\downarrow\\}}{\rm
U}{n}_{u}{n}_{l\sigma},$ (1) $\displaystyle H_{\rm B}$ $\displaystyle=$
$\displaystyle H_{\rm L}+H_{\rm M}+H_{\rm R},$ (2) $\displaystyle=$
$\displaystyle\sum_{\sigma,k}(\epsilon_{{\rm L}\sigma k}-\mu_{{\rm
L}\sigma})b^{\dagger}_{{\rm L}\sigma k}b_{{\rm L}\sigma k}$ $\displaystyle+$
$\displaystyle\sum_{k}(\epsilon_{{\rm M}k}-\mu_{\rm M})b^{\dagger}_{{\rm
M}k}b_{{\rm M}k}+\sum_{q}\epsilon_{{\rm R}q}a^{\dagger}_{{\rm R}q}a_{{\rm
R}q},$ $\displaystyle H_{\rm I}$ $\displaystyle=$ $\displaystyle H_{\rm
IL}+H_{\rm IM}+H_{\rm IR},$ (3) $\displaystyle=$
$\displaystyle\hbar\sum_{\sigma,k}(t_{{\rm L}k}b^{\dagger}_{{\rm L}\sigma
k}d_{l\sigma}+t^{*}_{{\rm L}k}d^{\dagger}_{l\sigma}b_{{\rm L}\sigma k})$
$\displaystyle+$ $\displaystyle\hbar\sum_{k}(t_{{\rm M}k}b^{\dagger}_{{\rm
M}k}d_{u}+t^{*}_{{\rm M}k}d^{\dagger}_{u}b_{{\rm M}k})$ $\displaystyle+$
$\displaystyle\hbar\sum_{q}(g_{{\rm R}q}a^{\dagger}_{{\rm
R}q}d^{\dagger}_{l\uparrow}d_{l\downarrow}+g^{*}_{{\rm
R}q}d^{\dagger}_{l\downarrow}d_{l\uparrow}a_{{\rm R}q}).$

Equation (1) represents the total system Hamiltonian of the two Coulomb-
coupled QDs, where ${\rm U}$ describes the long-range positive Coulomb
repulsion energy that permits energy exchange but forbids any particle
exchange between the QDs. The operator
${n}_{l\sigma}=d^{\dagger}_{l\sigma}d_{l\sigma}$ is the number operator for
${\rm QD}_{l}$, with eigenstates
$|\phi_{l}\rangle=\\{|0\rangle,\ket{\uparrow},\ket{\downarrow}\\}$ and
corresponding eigenenergies $0$, $\varepsilon_{l\uparrow}$ and
$\varepsilon_{l\downarrow}$, respectively, where $d^{\dagger}_{l\sigma}$
($d_{l\sigma}$) denotes the electron creation (annihilation) operator with a
single particle energy level $\varepsilon_{l\sigma}$, obeying anti-commutation
relation
$\\{d_{l\sigma},d^{\dagger}_{l\sigma^{\prime}}\\}=\delta_{\sigma\sigma^{\prime}}$;
$\sigma$ being the spin orientation of the electrons. Similarly,
$n_{u}=d^{\dagger}_{u}d_{u}$ is the number operator for ${\rm QD}_{u}$, with
eigenstates $|\phi_{u}\rangle={\\{|0\rangle,|1\rangle}\\}$ and corresponding
eigenenergies $0$ and $\varepsilon_{u}$, respectively, where,
$d^{\dagger}_{u}$ ($d_{u}$) represents the electron creation (annihilation)
operator for ${\rm QD}_{u}$, with a single particle energy level of
$\varepsilon_{u}$, satisfying the anti-commutation relation
$\\{d_{u},d^{\dagger}_{u}\\}=1$.

Equation (2) describes the total bath Hamiltonian $H_{\rm B}$, wherein $H_{\rm
L}$, $H_{\rm M}$ and $H_{\rm R}$ are the respective Hamiltonians for the left
(L), middle (M) and right (R) reservoirs. The operators $b^{\dagger}_{{\rm
L}\sigma k}$ ($b^{\dagger}_{{\rm M}k}$) and $b_{{\rm L}\sigma k}$ ($b_{{\rm
M}k}$) represent the creation and annihilation operators of electrons for the
L and M baths, where, $\epsilon_{{\rm L}\sigma k}$ and $\mu_{{\rm L}\sigma}$
stand for the energy and chemical potential of electrons corresponding to the
spinful fermionic reservoir (L), with $k$ being the continuous wave number
(momentum) and $\sigma$ denotes the electron spin. The difference between the
chemical potentials $\mu_{\rm L\downarrow}$ and $\mu_{\rm L\uparrow}$ is given
by the spin bias voltage i.e., $\Delta\mu_{\rm S}=\mu_{\rm
L\downarrow}-\mu_{\rm L\uparrow}$. On the other hand, $\epsilon_{{\rm M}k}$
and $\mu_{\rm M}$ refer to the energy and chemical potential of electrons
without any specific spin orientation for the spinless fermionic reservoir
(M). For the magnon bath (R), $a^{\dagger}_{{\rm R}q}$ and $a_{{\rm R}q}$ are
the bosonic creation and annihilation operators with the energy
$\epsilon_{{\rm R}q}$ and momentum $q$ respectively.

Equation (3) provides the total system-reservoir interaction Hamiltonian
$H_{\rm I}$, where $H_{\rm I\alpha}$ ($\alpha={\rm L,M,R}$) represents the
interaction between the system and the $\alpha$-th reservoir. Here the ${\rm
QD}_{l}$ (${\rm QD}_{u}$) is tunnel-coupled to the L and M reservoir with the
tunneling amplitudes $t_{\rm L(M)}$, allowing both particle and energy
exchange with the QDs, while ${\rm QD}_{l}$ is simultaneously coupled to a
magnon bath which flips only one spin at a time. Under strong coupling, the
eigenstates of $H_{\rm S}$ are determined by the tensor product of the number
operator’s eigenbasis $|\phi_{u}\phi_{l}\rangle$ of the coupled QD system. For
convenience, the six microstates of the coupled system
${\\{|0\rangle,|1\rangle}\\}\otimes\\{|0\rangle,\ket{\uparrow},\ket{\downarrow}\\}$,
are labeled by $|\mathbb{1}\rangle=|00\rangle$,
$|\mathbb{2}\rangle=|10\rangle$, $|\mathbb{3}\rangle=\ket{0\uparrow}$,
$|\mathbb{4}\rangle=\ket{0\downarrow}$, $|\mathbb{5}\rangle=\ket{1\uparrow}$,
$|\mathbb{6}\rangle=\ket{1\downarrow}$ and their corresponding eigenenergies
($\varepsilon_{\mathbb{i}}$, $\mathbb{i=1,2,....,6}$) are given by
$\varepsilon_{\mathbb{1}}=0$, $\varepsilon_{\mathbb{2}}=\varepsilon_{u}$,
$\varepsilon_{\mathbb{3}}=\varepsilon_{l\uparrow}$,
$\varepsilon_{\mathbb{4}}=\varepsilon_{l\downarrow}$,
$\varepsilon_{\mathbb{5}}=\varepsilon_{u}+\varepsilon_{l\uparrow}+{\rm U}$ and
$\varepsilon_{\mathbb{6}}=\varepsilon_{u}+\varepsilon_{l\downarrow}+{\rm U}$
respectively. There are in total nine allowed transitions: The transitions
$|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle$,
$|\mathbb{1}\rangle\leftrightarrow|\mathbb{4}\rangle$,
$|\mathbb{2}\rangle\leftrightarrow|\mathbb{5}\rangle$ and
$|\mathbb{2}\rangle\leftrightarrow|\mathbb{6}\rangle$ are driven by the
reservoir L, while the transitions
$|\mathbb{1}\rangle\leftrightarrow|\mathbb{2}\rangle$,
$|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle$ and
$|\mathbb{4}\rangle\leftrightarrow|\mathbb{6}\rangle$ are induced by the
reservoir M, and the transitions
$|\mathbb{3}\rangle\leftrightarrow|\mathbb{4}\rangle$ and
$|\mathbb{5}\rangle\leftrightarrow|\mathbb{6}\rangle$ are triggered by the
bath R.

To calculate the thermal spin ($J_{\rm S}$) and energy current ($J_{\rm E}$)
under the SSE and SPE, we first derive the Lindblad quantum master equation of
the reduced density matrix $\rho$ for the coupled QDs system under the Born-
Markov and Secular (BMS) approximation Breuer and Petruccione (2007);
Strasberg (2022); Gupt  _et al._ (2022, 2023) (see Appendix A)

$\frac{d\rho}{dt}=\mathcal{L}_{\rm L}[\rho]+\mathcal{L}_{\rm
R}[\rho]+\mathcal{L}_{\rm M}[\rho].$ (4)

Here $\mathcal{L}_{\alpha}$ ($\alpha={\rm L,R,M}$) is the Lindbladian due to
the interaction of the quantum system with its $\alpha$-th reservoir. The
explicit form of the superoperator $\mathcal{L}$ is given in terms of
dissipater

$\mathcal{D}(C)[\rho]=C\rho
C^{\dagger}-\frac{1}{2}\\{\rho,C^{\dagger}C\\},\;C\in\\{d_{l\sigma},d_{u},d^{\dagger}_{l\uparrow}d_{l\downarrow}\\},$
(5)

as follows:

$\displaystyle\mathcal{L}_{\rm L}[\rho]$ $\displaystyle=$
$\displaystyle\sum_{\sigma=\\{\uparrow,\downarrow\\}}\mathcal{L}_{{\rm
L}\sigma}[\rho],$ $\displaystyle\mathcal{L}_{{\rm L}\sigma}[\rho]$
$\displaystyle=$ $\displaystyle\sum_{\\{\varepsilon_{{\rm
L}\sigma}\\}}\gamma_{\rm L}\Big{[}f(\varepsilon_{{\rm L}\sigma},\mu_{{\rm
L}\sigma},T_{\rm L})\mathcal{D}(d^{\dagger}_{l\sigma})[\rho]$ (6)
$\displaystyle+$ $\displaystyle(1-f(\varepsilon_{{\rm L}\sigma},\mu_{{\rm
L}\sigma},T_{\rm L}))\mathcal{D}(d_{l\sigma})[\rho]\Big{]},$
$\displaystyle\mathcal{L}_{\rm M}[\rho]$ $\displaystyle=$
$\displaystyle\sum_{\\{\varepsilon_{\rm M}\\}}\gamma_{\rm
M}\Big{[}f(\varepsilon_{\rm M},\mu_{\rm M},T_{\rm
M})\mathcal{D}(d^{\dagger}_{u})[\rho]$ (7) $\displaystyle+$
$\displaystyle(1-f(\varepsilon_{\rm M},\mu_{\rm M},T_{\rm
M}))\mathcal{D}(d_{u})[\rho]\Big{]},$ $\displaystyle\mathcal{L}_{\rm R}[\rho]$
$\displaystyle=$ $\displaystyle\sum_{\\{\varepsilon_{\rm R}\\}}\gamma_{\rm
R}\Big{[}n(\varepsilon_{\rm R},T_{\rm
R})\mathcal{D}(d^{\dagger}_{l\downarrow}d_{l\uparrow})[\rho]$ (8)
$\displaystyle+$ $\displaystyle(1+n(\varepsilon_{\rm R},T_{\rm
R}))\mathcal{D}(d^{\dagger}_{l\uparrow}d_{l\downarrow})[\rho]\Big{]}.$

Note that we have implemented the strong coupling formalism to derive the
interaction picture master equation presented above Werlang  _et al._ (2014);
Ghosh  _et al._ (2022). Here, the strong coupling refers to the interaction
between the two QDs, while the system-reservoir coupling is assumed to be
weak, allowing for the safe implementation of the BMS approximation. In Eqs.
(6)-(8), all $\gamma$ values stand for the bare tunneling rates associated
with individual processes and depend on the system-reservoir coupling strength
through the respective bath spectral function. Lastly,
$f(\varepsilon,\mu,T)=[e^{(\varepsilon-\mu)/k_{B}T}+1]^{-1}$ and
$n(\varepsilon,T)=[e^{\varepsilon/k_{B}T}-1]^{-1}$ are respectively the Fermi-
Dirac (FD) and Bose-Einstein (BE) distribution functions with the positive
transition energy $\varepsilon$, chemical potential $\mu$ and temperature $T$
associated with the thermal reservoir, where $k_{B}$ is the Boltzmann
constant. Since the Hamiltonian $H_{\rm S}$ in Eq. (1) is diagonal in the
number state eigenbasis of the coupled QDs system, the reduced density matrix
$\rho$ of the above Lindblad master equation effectively decouples the
diagonal and off-diagonal matrix elements in the eigenbasis of $H_{\rm S}$
Ghosh  _et al._ (2022). The diagonal elements of the density matrix $\rho$
signify the occupation probabilities of each microstate and the time evolution
is given by

$\displaystyle\frac{dP_{\mathbb{1}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{12}}+J_{\mathbb{13}}+J_{\mathbb{14}},$ (9)
$\displaystyle\frac{dP_{\mathbb{2}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{21}}+J_{\mathbb{25}}+J_{\mathbb{26}},$ (10)
$\displaystyle\frac{dP_{\mathbb{3}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{31}}+J_{\mathbb{34}}+J_{\mathbb{35}},$ (11)
$\displaystyle\frac{dP_{\mathbb{4}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{41}}+J_{\mathbb{43}}+J_{\mathbb{46}},$ (12)
$\displaystyle\frac{dP_{\mathbb{5}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{52}}+J_{\mathbb{53}}+J_{\mathbb{56}},$ (13)
$\displaystyle\frac{dP_{\mathbb{6}}}{dt}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{62}}+J_{\mathbb{64}}+J_{\mathbb{65}}.$ (14)

Here $J_{\mathbb{ij}}$ stands for the net transition rate from state
$|\mathbb{j}\rangle$ to $|\mathbb{i}\rangle$ which is given by

$\displaystyle J_{\mathbb{ij}}$ $\displaystyle=$ $\displaystyle
k_{\mathbb{ij}}P_{\mathbb{j}}-k_{\mathbb{ji}}P_{\mathbb{i}},$ (15)
$\displaystyle J_{\mathbb{ij}}$ $\displaystyle=$ $\displaystyle-
J_{\mathbb{ji}},\;\quad\mathbb{i,j=1,2,....,6}$ (16)

where, $P_{\mathbb{i}}=\langle i|\rho|i\rangle$ is the population of the
$\mathbb{i}$-th eiegenstate and $k_{\mathbb{ji}}$
($k_{\mathbb{\ket{j}\leftarrow\ket{i}}}$) gives the transition probability
from microstate $|\mathbb{i}\rangle$ to microstate $|\mathbb{j}\rangle$. The
rate expressions $k_{\mathbb{ji}}$ for all transitions in terms of $\gamma$
and the distribution functions can be summarized as follows:

$\displaystyle k_{\mathbb{31}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}f(\varepsilon_{l\uparrow},\mu_{{\rm L}\uparrow},T_{\rm L}),$ $\displaystyle
k_{\mathbb{13}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}[1-f(\varepsilon_{l\uparrow},\mu_{{\rm L}\uparrow},T_{\rm L})],$
$\displaystyle k_{\mathbb{41}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}f(\varepsilon_{l\downarrow},\mu_{{\rm L}\downarrow},T_{\rm L}),$
$\displaystyle k_{\mathbb{14}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}[1-f(\varepsilon_{l\downarrow},\mu_{{\rm L}\downarrow},T_{\rm L})],$
$\displaystyle k_{\mathbb{52}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}f(\varepsilon_{l\uparrow}+{\rm U},\mu_{{\rm L}\uparrow},T_{\rm L}),$
$\displaystyle k_{\mathbb{25}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}[1-f(\varepsilon_{l\uparrow}+{\rm U},\mu_{{\rm L}\uparrow},T_{\rm L})],$
$\displaystyle k_{\mathbb{62}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}f(\varepsilon_{l\downarrow}+{\rm U},\mu_{{\rm L}\downarrow},T_{\rm L}),$
$\displaystyle k_{\mathbb{26}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
L}[1-f(\varepsilon_{l\downarrow}+{\rm U},\mu_{{\rm L}\downarrow},T_{\rm L})],$
$\displaystyle k_{\mathbb{21}}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
M}f(\varepsilon_{u},\mu_{\rm M},T_{\rm M}),$ $\displaystyle k_{\mathbb{12}}$
$\displaystyle=$ $\displaystyle\gamma_{\rm M}[1-f(\varepsilon_{u},\mu_{\rm
M},T_{\rm M})],$ $\displaystyle k_{\mathbb{53}}$ $\displaystyle=$
$\displaystyle k_{\mathbb{64}}=\gamma_{\rm M}f(\varepsilon_{u}+{\rm
U},\mu_{\rm M},T_{\rm M}),$ $\displaystyle k_{\mathbb{35}}$ $\displaystyle=$
$\displaystyle k_{\mathbb{46}}=\gamma_{\rm M}[1-f(\varepsilon_{u}+{\rm
U},\mu_{\rm M},T_{\rm M})],$ $\displaystyle k_{\mathbb{43}}$ $\displaystyle=$
$\displaystyle k_{\mathbb{65}}=\gamma_{\rm
R}n(\varepsilon_{l\downarrow}-\varepsilon_{l\uparrow},T_{\rm R}),$
$\displaystyle k_{\mathbb{34}}$ $\displaystyle=$ $\displaystyle
k_{\mathbb{56}}=\gamma_{\rm
R}[1+n(\varepsilon_{l\downarrow}-\varepsilon_{l\uparrow},T_{\rm R})].$ (17)

Combining Eqs. (9)-(14) with Eq. (15), it is evident that the evolution
equations for the microscopic probabilities exhibit linearity with respect to
the populations $\\{P_{\mathbb{i}}\\}$. As a result, we can cast these
equations in the following compact form

$\frac{dP_{\mathbb{i}}}{dt}=\sum^{\mathbb{6}}_{\mathbb{j=1}}J_{\mathbb{ij}}=\sum^{\mathbb{6}}_{\mathbb{j=1}}k_{\mathbb{ij}}P_{\mathbb{j}}-k_{\mathbb{ji}}P_{\mathbb{i}};\quad\mathbb{i\neq
j},$ (18)

where $\sum^{\mathbb{6}}_{\mathbb{i=1}}P_{\mathbb{i}}=1$. Equation (18) is
known as the quantum kinetic Pauli master equation which is “classical” in
looking but quantum mechanical in content through the transition probabilities
$\\{k_{\mathbb{ij}}\\}$, determined by the Fermi’s golden rule within BMS
approximation and the statistical properties of the respective quantum baths
Sinha  _et al._ (2011a, b). To obtain the steady-state solution
$\bar{P}_{\mathbb{i}}$ of Eq.(18), one has to solve the system of linear
equations, satisfying the conditions $0\leq\bar{P}_{\mathbb{i}}\leq 1$ and
$\sum_{\mathbb{i}}\bar{P}_{\mathbb{i}}=1$. With the help of Eq. (18) and the
rate coefficients calculated from the above microscopic picture (Cf. (II)), it
is possible to evaluate the steady-state spin and energy currents in terms of
the net transition rates (Cf. Eq. (15)) between the system microstates, where
$\\{P_{\mathbb{i}}\\}$ get replaced by the steady-state populations
$\\{\bar{P}_{\mathbb{i}}\\}$. Following the definition of the spin and energy
currents mentioned in Appendix-A, we obtain the mathematical expression of the
steady-state spin current $J_{\rm S}$ which flows from left to right, as

$\displaystyle J_{\rm S}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}\Tr\\{d^{\dagger}_{l\downarrow}d_{l\downarrow}\mathcal{L}_{\rm
L\downarrow}[\rho]\\}-\Tr\\{d^{\dagger}_{l\uparrow}d_{l\uparrow}\mathcal{L}_{\rm
L\uparrow}[\rho]\\}\Big{)},$ $\displaystyle=$ $\displaystyle
J_{\mathbb{34}}+J_{\mathbb{56}}=(k_{\mathbb{34}}\bar{P}_{\mathbb{4}}-k_{\mathbb{43}}\bar{P}_{\mathbb{3}})+(k_{\mathbb{56}}\bar{P}_{\mathbb{6}}-k_{\mathbb{65}}\bar{P}_{\mathbb{5}}),$

and the steady-state energy (heat) current $J_{\rm E}$, through the middle
reservoir is given by

$\displaystyle J_{\rm E}$ $\displaystyle=$
$\displaystyle\Tr\\{\mathcal{L}_{\rm M}[\rho]H_{\rm S}\\}={\rm
U}(J_{\mathbb{53}}+J_{\mathbb{64}})$ (20) $\displaystyle=$ $\displaystyle{\rm
U}[(k_{\mathbb{53}}\bar{P}_{\mathbb{3}}-k_{\mathbb{35}}\bar{P}_{\mathbb{5}})+(k_{\mathbb{64}}\bar{P}_{\mathbb{4}}-k_{\mathbb{46}}\bar{P}_{\mathbb{6}})].$

Eq. (20) immediately implies that a finite energy current always requires a
finite Coulomb interaction energy. However, obtaining the exact analytical
solutions for $J_{\rm S}$, $J_{\rm E}$ in terms of steady-state populations
[Eqs. (II) and (20)] by solving the linear master equation [Eq. (18)] is by no
means a trivial task. Secondly, while, one may in principle use exact Eqs.
(LABEL:spin-current-Js) and (20) to numerically compute the steady-state spin
and energy currents, it does not provide any physical insight into the
underlying transport mechanisms leading to SSE and SPE. Nor does it explain
how the macroscopic spin and energy currents are related to the thermodynamic
forces that give rise to spin-thermoelectric effects as a manifestation of
thermodynamic cross-effects.

An alternative yet effective method is to calculate algebraic expressions for
steady-state currents through a network or mathematical graph theory Tutte
(2001); Balakrishnan and Ranganathan (2012). This also allows us to understand
the operational principles of QD-based spin-thermoelectric effects quite
easily. In this method, one first constructs a basic graph $\mathbb{G}$ as a
diagrammatic representation of the right-hand side of Eq. (18). To extract the
principal mechanism from complex transport behaviors, one then decomposes the
quantum transition network into cycle trajectories, collects the cycle fluxes
using algebraic graph theory, and selects the top-ranked cycle fluxes—i.e.,
the cycle trajectories with the highest probabilities Wang  _et al._ (2022).
In the following section, we illustrate this method in the context of the
present problem and establish the connection between the microscopic
descriptions of the non-equilibrium system via the basic graph and the
macroscopic description of thermolectric phenomena in terms of thermodynamic
forces and fluxes, including the celebrated Onsager and Kelvin relations
Callen (1985). The key concept throughout the entire formalism is the
expression of entropy production rate in the framework of graph theory Landi
and Paternostro (2021); Schnakenberg (1976).

##  III Network theory and Reciprocity relation

The network or graph theory found its first application in electricity, with
Kirchhoff making a pioneering contribution to the understanding of electrical
circuits as non-equilibrium systems involving electric current and potential.
Since then, graph theory has expanded its horizons and produced a flurry of
inspiring early works by Hill, Kohler, Vollmerhaus, King, and Altman Hill and
Chen (1975); Kohler and Vollmerhaus (1980); King and Altman (1956),
particularly on biophysical and biochemical systems. A vast body of literature
is available on this subject Wu  _et al._ (2012); Einax  _et al._ (2011);
Einax and Nitzan (2014); Ren (2017); Dutta  _et al._ (2020); still,
Schnakenberg’s 1976 review is considered a seminal contribution to this field
Schnakenberg (1976).

###  III.1 Quantum Transition network and Cycle flux analysis

As an extension of network theory to quantum systems, the notable work of Wang
et al. Wang  _et al._ (2022) is worth mentioning. They have recently
demonstrated that the dissipative quantum dynamics of non-equilibrium
transport can be mapped onto networks of quantum state transitions, where
nodes or vertices correspond to quantum states, and the connecting lines or
edges between two quantum states represent their allowed transitions.

Figure 2: Schematic diagram of the basic graph ($\mathbb{G}$). Subcycles
$\\{\mathcal{C}_{1},\mathcal{C}_{5},\mathcal{C}_{8},\mathcal{C}_{9},\mathcal{C}_{10}\\}$,
sharing the common edge ($\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{4}}$)
are used to calculate edge flux $J_{\mathbb{43}}$ [Cf. Eq.(24)].

In the present case, the diagrammatic representation of the quantum transport
processes under the non-equilibrium condition is shown in Fig. 2 in the form
of a basic graph ($\mathbb{G}$), where each node or vertex represents a
quantum state $\\{|\mathbb{i}\rangle\\}$ along with its associated
(microscopic) occupation probability $\\{P_{\mathbb{i}}\\}$. The transition
between adjacent quantum states $|\mathbb{i}\rangle$ and $|\mathbb{j}\rangle$
are depicted by edges. The steady-state population $\bar{P}_{\mathbb{i}}$ can
then be calculated as

$\bar{P}_{\mathbb{i}}=\frac{\Lambda_{\mathbb{i}}}{\Lambda},\quad\text{with}\quad
0\leq\bar{P}_{\mathbb{i}}\leq
1\quad\text{and}\quad\sum_{\mathbb{i}}\bar{P}_{\mathbb{i}}=1;$ (21)

where $\Lambda_{\mathbb{i}}$ represents the sum of the weight of the spanning
trees rooted on the $\mathbb{\ket{i}}$-th state and $\Lambda$ is defined as
the sum of the weights of the spanning trees rooted on every individual state
$\\{\mathbb{\ket{i}}\\}$, i.e. $\sum_{\mathbb{i}}\Lambda_{\mathbb{i}}$. In the
literature, the above method is known as Kirchhoff’s theorem Schnakenberg
(1976); Kirchhoff (1847). According to this theorem, a spanning tree is a
subgraph of $\mathbb{G}$ which includes all the vertices with the minimum
number of edges that are always connected but have no circuits (cyclic
sequence of edges or cycle trajectory). To construct a spanning tree, one
should remove $\nu=e-v+1$ number of edges of the basic graph $\mathbb{G}$,
where $e$ and $v$ are the numbers of edges and vertices in $\mathbb{G}$
Schnakenberg (1976). As a result, all possible spanning trees contain an equal
number of vertices and edges.

Under the non-equilibrium condition, each edge represents a transport process
and the rate of these transport processes is determined by the net transition
rate or edge flux. The steady-state edge flux from a state $\mathbb{\ket{j}}$
and $\mathbb{\ket{i}}$ is defined as

$J_{\mathbb{ij}}=k_{\mathbb{ij}}\bar{P}_{\mathbb{j}}-k_{\mathbb{ji}}\bar{P}_{\mathbb{i}},$
(22)

where, each edge denotes a pair of transitions with the transition
probabilities $k_{\mathbb{ij}}$ (from $|\mathbb{j}\rangle$ to
$|\mathbb{i}\rangle$) and $k_{\mathbb{ji}}$ (from $|\mathbb{i}\rangle$ to
$|\mathbb{j}\rangle$) Schnakenberg (1976); Wu  _et al._ (2012). Typically, a
basic graph $\mathbb{G}$ is comprised of numerous undirected subcycles
($\mathcal{C}$), and each of these subcycles represents a pair of two one-
directional circuits [Fig. 2], namely $\mathcal{C}^{+}$ (counterclockwise) and
$\mathcal{C}^{-}$ (clockwise) Kohler and Vollmerhaus (1980). Since the
circuits are formed by the cyclic sequence of edges within $\mathbb{G}$, the
edge flux can be defined in terms of the circuit fluxes Schnakenberg (1976),
as

$J_{\mathbb{ij}}=\sum_{\mathcal{C}}\mathcal{S}_{ij}(\mathcal{C})(J^{+}_{\mathcal{C}}-J^{-}_{\mathcal{C}})=\sum_{\mathcal{C}}\mathcal{S}_{ij}(\mathcal{C})J_{\mathcal{C}}.$
(23)

Here $J_{\mathcal{C}}=J^{+}_{\mathcal{C}}-J^{-}_{\mathcal{C}}$ denotes the net
cycle flux wherein $J^{+}_{\mathcal{C}}$ and $J^{-}_{\mathcal{C}}$ are the
circuit fluxes correspond to circuits $\mathcal{C}^{+}$ and $\mathcal{C}^{-}$
respectively, with the prefactor $\mathcal{S}_{\mathbb{ij}}(\mathcal{C})=0,\pm
1$. $\mathcal{S}_{\mathbb{ij}}(\mathcal{C})=0$ if $\mathcal{C}^{+}$ and
$\mathcal{C}^{-}$ does not contain the edge
$\mathbb{\ket{j}}\rightarrow\mathbb{\ket{i}}$;
$\mathcal{S}_{\mathbb{ij}}(\mathcal{C})=+1$ if the orientation of
$\mathcal{C}^{+}$ ($\mathcal{C}^{-}$) is along (opposite) to edge
$\mathbb{\ket{j}}\rightarrow\mathbb{\ket{i}}$ and
$\mathcal{S}_{\mathbb{ij}}(\mathcal{C})=-1$ if the orientation of
$\mathcal{C}^{+}$ ($\mathcal{C}^{-}$) is opposite (along) to edge
$\mathbb{\ket{j}}\rightarrow\mathbb{\ket{i}}$. For example, the edge flux
$J_{\mathbb{43}}$ ($J_{\mathbb{\ket{4}\leftarrow\ket{3}}}$) in the basic graph
$\mathbb{G}$, can be expressed in terms of the circuit fluxes [Fig. 2] as

$J_{\mathbb{43}}=J^{+}_{\mathcal{C}_{1}}-J^{-}_{\mathcal{C}_{1}}-J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{5}}+J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}}-J^{+}_{\mathcal{C}_{9}}+J^{-}_{\mathcal{C}_{9}}-J^{+}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{10}}.$
(24)

The name “circuit” was initially introduced by Kohler and Vollmerhaus Kohler
and Vollmerhaus (1980) and also termed a “one-way cycle” by Hill Hill and
Kedem (1966). However, we prefer to use the term “circuit” or “cycle
trajectory” to avoid confusion with the usual “cycle”. Hill and Chen provided
the physical interpretation for circuit fluxes Hill and Chen (1975), revealing
that these fluxes signify the ‘frequency’ (or rate) of circuit completions
along a particular cycle trajectory. To be specific, the circuit flux
associated with a one-directional cycle trajectory $\mathcal{C^{\pm}}$ is
given by

$J^{\pm}_{\mathcal{C}}=\Pi^{\pm}_{\mathcal{C}}\frac{\Lambda_{\mathcal{C}}}{\Lambda}.$
(25)

Here, $\Pi^{\pm}_{\mathcal{C}}$ denotes the weight factor which is determined
by the product of the transition rates along the circuit $\mathcal{C^{\pm}}$.
For example, the clockwise cycle trajectory
$\mathcal{C}^{-}_{1}(\mathbb{\ket{1}}\rightarrow\mathbb{\ket{4}}\rightarrow\mathbb{\ket{3}}\rightarrow\mathbb{\ket{1}})$
[Fig. 2] has the weight factor
$\Pi^{-}_{\mathcal{C}_{1}}=k_{\mathbb{13}}k_{\mathbb{34}}k_{\mathbb{41}}$,
where, $\Lambda_{\mathcal{C}}$ represents the sum of the weight of the
spanning trees rooted on cycle $\mathcal{C}$ and
$\Lambda=\sum_{\mathbb{i}}\Lambda_{\mathbb{i}}$.

Figure 3: Spanning trees rooted on cycle $\mathcal{C}_{3}$ (shaded region) of
the basic graph.

Now, there are a total of 22 paired cycle trajectories, or 11 subcycles, for
our basic graph $\mathbb{G}$, which are as follows:

$\displaystyle\mathcal{C}_{1}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{2}$ $\displaystyle:$
$\displaystyle|\mathbb{2}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{2}\rangle$
$\displaystyle\mathcal{C}_{3}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{4}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{5}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{6}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{7}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{8}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{1}\rangle$
$\displaystyle\mathcal{C}_{9}$ $\displaystyle:$
$\displaystyle|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{3}\rangle$
$\displaystyle\mathcal{C}_{10}$ $\displaystyle:$
$\displaystyle|\mathbb{2}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{2}\rangle$
$\displaystyle\mathcal{C}_{11}$ $\displaystyle:$
$\displaystyle|\mathbb{1}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{6}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{1}\rangle.$
(26)

Hence, enumerating a large number of spanning trees rooted at each individual
state, as well as for cycles, poses a formidable challenge. This difficulty
becomes more pronounced with the increasing size of the basic graph.

To bypass this problem, we utilize the generalized matrix-tree theorem from
algebraic graph Wang  _et al._ (2022); Gupt  _et al._ (2023) by rewriting the
master equation in the following form Keizer (1972):
$\dot{\rm\textbf{P}}=-{\rm\textbf{MP}}$, where
${\rm\textbf{P}}=\\{P_{\mathbb{1}},P_{\mathbb{2}},P_{\mathbb{3}},P_{\mathbb{4}},P_{\mathbb{5}},P_{\mathbb{6}}\\}$
is a column matrix and M is a square matrix, given by

$\displaystyle{\rm\textbf{M}}=\left[{\begin{array}[]{cccccc}k_{\mathbb{21}}+k_{\mathbb{31}}+k_{\mathbb{41}}&-k_{\mathbb{21}}&k_{\mathbb{13}}&-k_{\mathbb{14}}&0&0\\\
-k_{\mathbb{21}}&k_{\mathbb{12}}+k_{\mathbb{52}}+k_{\mathbb{62}}&0&0&-k_{\mathbb{25}}&-k_{\mathbb{26}}\\\
-k_{\mathbb{31}}&0&k_{\mathbb{13}}+k_{\mathbb{43}}+k_{\mathbb{53}}&-k_{\mathbb{34}}&-k_{\mathbb{35}}&0\\\
-k_{\mathbb{41}}&0&-k_{\mathbb{43}}&k_{\mathbb{14}}+k_{\mathbb{34}}+k_{\mathbb{64}}&0&-k_{\mathbb{46}}\\\
0&-k_{\mathbb{52}}&-k_{\mathbb{53}}&0&k_{\mathbb{25}}+k_{\mathbb{35}}+k_{\mathbb{65}}&-k_{\mathbb{56}}\\\
0&-k_{\mathbb{62}}&0&-k_{\mathbb{64}}&-k_{\mathbb{65}}&k_{\mathbb{26}}+k_{\mathbb{46}}+k_{\mathbb{56}}\\\
\end{array}}\right].$ (33)

Equation (33) is known as the Laplacian or transition matrix of the weighted
graph $\mathbb{G}$. Furthermore, in accordance with the matrix tree theorem,
it is possible to compute both the numerator and denominator of Eqs. (21) and
(25) as the determinants of the reduced transition matrix. For instance,
$\Lambda_{\mathbb{i}}$ is related to ${\rm\textbf{M}}[\mathbb{i},\mathbb{i}]$
which can be obtained by removing the $\mathbb{i}$-th row and column of the
Laplacian matrix M. Similarly, $\Lambda_{\mathcal{C}}$ is identical to the
$\det({\rm\textbf{M}}[\mathcal{C},\mathcal{C}])$, obtained by deleting rows
and columns belonging to cycle $\mathcal{C}$ of the transition matrix M. This
directly leads to a simple algebraic expression of the steady-state population

$\bar{P}_{\mathbb{i}}=\frac{\det({\rm\textbf{M}}[\mathbb{i};\mathbb{i}])}{\sum_{\mathbb{i}}\det({\rm\textbf{M}}[\mathbb{i};\mathbb{i}])},$
(34)

and the one-directional circuit flux associated with circuits
$\mathcal{C}^{\pm}$ in the following form

$J^{\pm}_{\mathcal{C}}=\Pi^{\pm}_{\mathcal{C}}\frac{\det({\rm\textbf{M}}[\mathcal{C};\mathcal{C}])}{\sum_{\mathbb{i}}\det({\rm\textbf{M}}[\mathbb{i};\mathbb{i}])},$
(35)

where
${\sum_{\mathbb{i}}\det({\rm\textbf{M}}[\mathbb{i};\mathbb{i}])}=\Lambda$. As
an example, the sum of the weights of spanning trees rooted on cycle
$\mathcal{C}_{3}$ in terms of the reduced determinant of the original
Laplacian matrix M is obtained by deleting rows and columns
$\mathbb{i(1,2,3,5)}\in\mathcal{C}_{3}$ [Fig. 3]. So, the principle minor
${\rm\textbf{M}}[\mathcal{C}_{3},\mathcal{C}_{3}]$ or
${\rm\textbf{M}}[1,3,5,2;1,3,5,2]$ and its determinant takes the form of

$\displaystyle{\rm\textbf{M}}[\mathcal{C}_{3},\mathcal{C}_{3}]=\left[{\begin{array}[]{ccc}k_{\mathbb{14}}+k_{\mathbb{34}}+k_{\mathbb{64}}&-k_{\mathbb{46}}\\\
-k_{\mathbb{64}}&k_{\mathbb{26}}+k_{\mathbb{46}}+k_{\mathbb{56}}\end{array}}\right],$
(38)

and

$\displaystyle\det({\rm\textbf{M}}[\mathcal{C}_{3},\mathcal{C}_{3}])=k_{\mathbb{14}}k_{\mathbb{26}}+k_{\mathbb{14}}k_{\mathbb{46}}+k_{\mathbb{14}}k_{\mathbb{56}}+k_{\mathbb{34}}k_{\mathbb{26}}$
$\displaystyle+k_{\mathbb{34}}k_{\mathbb{46}}+k_{\mathbb{34}}k_{\mathbb{56}}+k_{\mathbb{26}}k_{\mathbb{64}}+k_{\mathbb{56}}k_{\mathbb{64}}.$
(39)

respectively. Equation (III.1) contains a sum of eight terms, indicating that
there are eight possible spanning trees rooted on cycle $\mathcal{C}_{3}$.
Each term represents the weight of a spanning tree in which all the weighted
edges are directed towards the cycle $\mathcal{C}_{3}$ [Fig. 3]. Finally, from
Eq. (35), it is clear that we can calculate the circuit fluxes
$J^{\pm}_{\mathcal{C}}$ for each stochastic cycle trajectory with the
orientation either clockwise or counterclockwise and efficiently rank out the
top-ranked circuit fluxes. Next, we will show how the microscopic details of
cycle and circuit fluxes help us understand the SSE and SPE, connecting spin
and energy currents to macroscopic thermodynamic forces in the
phenomenological laws of irreversible thermodynamics.

###  III.2 Onsager Relation

The first step on our way from a microscopic to a macroscopic description is
to establish an expression for entropy production rate, the key quantity in
understanding any irreversible processes Schnakenberg (1976); Landi and
Paternostro (2021). To start with, we consider the von Neumann entropy

$\mathcal{S}=-k_{B}\sum_{\mathbb{i}}P_{\mathbb{i}}\ln P_{\mathbb{i}},$ (40)

in the framework of our discrete-state quantum transition network,
characterized by its microscopic probability $\\{P_{\mathbb{i}}\\}$. Thus, the
time evolution of $\mathcal{S}$ is given by

$\frac{d\mathcal{S}}{dt}=-k_{B}\sum_{\mathbb{i}}\frac{dP_{\mathbb{i}}}{dt}\ln
P_{\mathbb{i}}.$ (41)

With the help of the quantum kinetic Pauli master equation (18), one can
rewrite Eq. (41) as

$\displaystyle\frac{d\mathcal{S}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{2}k_{B}\sum_{\mathbb{i,j}}J_{\mathbb{ij}}\ln\Big(\frac{P_{\mathbb{j}}}{P_{\mathbb{i}}}\Big{missing}).$
(42)

Now, following Schnakenberg’s suggestion Schnakenberg (1976), we split Eq.
(42) into two parts,

$\frac{d\mathcal{S}}{dt}=\dot{\Phi}(t)+\dot{\sigma}(t),$ (43)

where we identify the first term $\dot{\Phi}(t)$ as the entropy flux rate

$\dot{\Phi}(t)=-\frac{1}{2}k_{B}\sum_{\mathbb{i,j}}J_{\mathbb{ij}}\ln\Big(\frac{k_{\mathbb{ij}}}{k_{\mathbb{ji}}}\Big{missing}),$
(44)

which arises from the interaction between the system and its surroundings. The
second term $\dot{\sigma}(t)$ is the total entropy production rate

$\displaystyle\dot{\sigma}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}k_{B}\sum_{\mathbb{i,j}}(k_{\mathbb{ij}}P_{\mathbb{j}}-k_{\mathbb{ji}}P_{\mathbb{i}})\ln\Big(\frac{k_{\mathbb{ij}}P_{\mathbb{j}}}{k_{\mathbb{ji}}P_{\mathbb{i}}}\Big{missing}).$
(45)

Equations (45) may appear a little artificial at first glance, and a natural
question to be raised at this point is whether Eq. (45) has anything to do
with the entropy production of the phenomenological irreversible
thermodynamics, which needs be expressed as a bilinear form of the macroscopic
thermodynamic forces and fluxes. It must be emphasized that neither
$\dot{\mathcal{S}}$ in Eq. (43) nor $\dot{\Phi}$ in Eq. (44) are necessarily
positive, but only $\dot{\sigma}(t)\geq 0$, since it takes a form
$(a-b)\ln(a/b)\geq 0$ [Cf. Eq. (45)]. This is indeed true since the total
entropy production rate ($\dot{\sigma}$) of any system must be always positive
Schnakenberg (1976); Landi and Paternostro (2021). Thus, it turns out that Eq.
(45) satisfies the basic criteria for the entropy production rate. Under the
steady-state condition, there is no change in the entropy of the system which
implies Schnakenberg (1976); Landi and Paternostro (2021)

$\dot{\sigma}=-\dot{\Phi}(t)=\frac{1}{2}k_{B}\sum_{\mathbb{i,j}}(k_{\mathbb{ij}}\bar{P}_{\mathbb{j}}-k_{\mathbb{ji}}\bar{P}_{\mathbb{i}})\ln\Big(\frac{k_{\mathbb{ij}}}{k_{\mathbb{ji}}}\Big{missing}).$
(46)

Using Eqs. (22) and (23), one may rewrite the above equation in terms of the
circuit and cycle fluxes as Schnakenberg (1976),

$\displaystyle\dot{\sigma}$ $\displaystyle=$ $\displaystyle
k_{B}\sum_{\mathcal{C}}(J^{+}_{\mathcal{C}}\mathcal{A}^{+}_{\mathcal{C}}+J^{-}_{\mathcal{C}}\mathcal{A}^{-}_{\mathcal{C}})$
(47) $\displaystyle=$ $\displaystyle
k_{B}\sum_{\mathcal{C}}(J^{+}_{\mathcal{C}}-J^{-}_{\mathcal{C}})\mathcal{A}^{+}_{\mathcal{C}}=k_{B}\sum_{\mathcal{C}}J_{\mathcal{C}}\mathcal{X}_{\mathcal{C}},$

where
$\mathcal{X}_{\mathcal{C}}=\mathcal{A}^{+}_{\mathcal{C}}=\ln(\Pi^{+}_{\mathcal{C}}/\Pi^{-}_{\mathcal{C}})=-\mathcal{A}^{-}_{\mathcal{C}}$
is called the the cycle affinity. For a given cycle $\mathcal{C}$, it measures
the imbalance or asymmetry between the transition rates along two opposite
cycle trajectories $\mathcal{C}^{\pm}$ and hence qualifies as a thermodynamic
force Ohga  _et al._ (2023). This is because, when
$\mathcal{X}_{\mathcal{C}}=0$, it implies $J_{\mathcal{C}}=0$, resulting in
equal circuit fluxes in both directions, i.e.,
$J^{+}_{\mathcal{C}}=J^{-}_{\mathcal{C}}$. Equation (47), expressed in terms
of cycle fluxes and cycle forces, can thus be regarded as a microscopic or
stochastic version of the phenomenological Onsager relation Landi and
Paternostro (2021).

Moreover, we find from Eq. (35) that the ratio of $J^{\pm}_{\mathcal{C}}$ is
equal to the ratio of weight factors $\Pi^{\pm}_{\mathcal{C}}$ for each cycle,
which, in turn, is determined by the ratio of the product of the transitions
rates along circuits $\mathcal{C}^{\pm}$ and can be computed in terms of
externally controllable, macroscopic physical quantities $T_{0}$, $\delta T$
and $\Delta\mu_{\rm S}$, as follows (see Appendix B):

$\displaystyle\frac{J^{+}_{\mathcal{C}_{1}}}{J^{-}_{\mathcal{C}_{1}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (48)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{2}}}{J^{-}_{\mathcal{C}_{2}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{2}}}{\Pi^{-}_{\mathcal{C}_{2}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (49)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{3}}}{J^{-}_{\mathcal{C}_{3}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{3}}}{\Pi^{-}_{\mathcal{C}_{3}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}},$ (50)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{4}}}{J^{-}_{\mathcal{C}_{4}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{4}}}{\Pi^{-}_{\mathcal{C}_{4}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (51)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{5}}}{J^{-}_{\mathcal{C}_{5}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{5}}}{\Pi^{-}_{\mathcal{C}_{5}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (52)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{6}}}{J^{-}_{\mathcal{C}_{6}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{6}}}{\Pi^{-}_{\mathcal{C}_{6}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}},$ (53)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{7}}}{J^{-}_{\mathcal{C}_{7}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{7}}}{\Pi^{-}_{\mathcal{C}_{7}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (54)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{8}}}{J^{-}_{\mathcal{C}_{8}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{8}}}{\Pi^{-}_{\mathcal{C}_{8}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (55)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{9}}}{J^{-}_{\mathcal{C}_{9}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{9}}}{\Pi^{-}_{\mathcal{C}_{9}}}=1,$
(56) $\displaystyle\frac{J^{+}_{\mathcal{C}_{10}}}{J^{-}_{\mathcal{C}_{10}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{10}}}{\Pi^{-}_{\mathcal{C}_{10}}}=e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (57)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{11}}}{J^{-}_{\mathcal{C}_{11}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{11}}}{\Pi^{-}_{\mathcal{C}_{11}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}},$ (58)

In order to derive Eqs. (48)-(58), one makes use of Eq. (II), where the
explicit form of the distribution functions are governed by the quantum
statistical properties of the respective thermal reservoirs (see Appendix B
for details). The advantage to writing the above set of equations as a ratio
of the circuits fluxes lies in its ability to indicate that if the external
biases $\Delta\mu_{\rm S}$ and $\delta T$ are zero on the r.h.s of Eqs.
(48)-(58), then regardless of the magnitudes the circuit fluxes, corresponding
cycle can’t contribute the spin and energy currents. Therefore, one can infer
the cycle fluxes associated with the subcycles
$\\{\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{10},\mathcal{C}_{11}\\}$ are
controlled by the spin bias voltage $\Delta\mu_{\rm S}$ and hence can only
contribute to the spin current $J_{\rm S}$. Whereas the net cycle fluxes
associated with the subcycles $\mathcal{C}_{3}$ and $\mathcal{C}_{6}$ are
dependent on the Coulomb interaction ${\rm U}$ and the temperature gradient
$\delta T$ and thereby contributing to the energy current $J_{\rm E}$.
However, there are few subcycles
$\\{\mathcal{C}_{4},\mathcal{C}_{5},\mathcal{C}_{7},\mathcal{C}_{8}\\}$ and
their conjugate fluxes are governed by both $\delta T$ as well as
$\Delta\mu_{\rm S}$ and therefore can contribute to both $J_{\rm S}$ and
$J_{\rm E}$. Indeed, these are the four cycles that are responsible for the
spin-thermoelectric cross-effects of SSE and SPE as we will demonstrate in
Sec. IV. Note that the net cycle flux associated with $\mathcal{C}_{9}$ is
identically zero because the circuit fluxes in both directions (clockwise and
counterclockwise) are equal as evident from Eq.(56). As a result, we can use
Eq. (23) to rewrite macroscopic spin and energy current expressions [Eqs. (II)
and (20)] in terms of the microscopic circuit and cycle fluxes as follows:

$\displaystyle J_{\rm S}$ $\displaystyle=$
$\displaystyle-(J^{+}_{\mathcal{C}_{1}}-J^{-}_{\mathcal{C}_{1}})-(J^{+}_{\mathcal{C}_{2}}-J^{-}_{\mathcal{C}_{2}})-(J^{+}_{\mathcal{C}_{4}}-J^{-}_{\mathcal{C}_{4}})$
$\displaystyle+$
$\displaystyle(J^{+}_{\mathcal{C}_{5}}-J^{-}_{\mathcal{C}_{5}})+(J^{+}_{\mathcal{C}_{7}}-J^{-}_{\mathcal{C}_{7}})-(J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}})$
$\displaystyle+$
$\displaystyle(J^{+}_{\mathcal{C}_{10}}-J^{-}_{\mathcal{C}_{10}})-(J^{+}_{\mathcal{C}_{11}}-J^{-}_{\mathcal{C}_{11}})$
$\displaystyle=$ $\displaystyle-
J_{\mathcal{C}_{1}}-J_{\mathcal{C}_{2}}-J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{7}}-J_{\mathcal{C}_{8}}$
$\displaystyle+$ $\displaystyle J_{\mathcal{C}_{10}}-J_{\mathcal{C}_{11}},$
(60) $\displaystyle J_{\rm E}$ $\displaystyle=$ $\displaystyle{\rm
U}[(J^{+}_{\mathcal{C}_{3}}-J^{-}_{\mathcal{C}_{3}})+(J^{+}_{\mathcal{C}_{4}}-J^{-}_{\mathcal{C}_{4}})+(J^{+}_{\mathcal{C}_{5}}-J^{-}_{\mathcal{C}_{5}})$
(61) $\displaystyle+$
$\displaystyle(J^{+}_{\mathcal{C}_{6}}-J^{-}_{\mathcal{C}_{6}})+(J^{+}_{\mathcal{C}_{7}}-J^{-}_{\mathcal{C}_{7}})+(J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}})]$
$\displaystyle=$ $\displaystyle{\rm
U}[J_{\mathcal{C}_{3}}+J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{6}}+J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}].$
(62)

Similar to Eqs. (II) and (20), the above set of equations are the most general
ones, however, the latter has an advantage over the previous set of equations.
Eqs. (61) and (62) can be expressed in terms of macroscopic forces,
facilitating the connection between SSE and SPE as a manifestation of
thermodynamic cross-effects. In order to identify the phenomenological forces,
we substitute Eqs. (48)-(58) into Eq. (47), to write the entropy production
rate as a sum over the associated thermodynamic forces and fluxes as (Appendix
B)

$\displaystyle\dot{\sigma}$ $\displaystyle=$ $\displaystyle{\rm
U}[J_{\mathcal{C}_{3}}+J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{6}}+J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}]\frac{\delta
T}{T_{0}(T_{0}+\delta T)}$ (63) $\displaystyle-$
$\displaystyle[J_{\mathcal{C}_{1}}+J_{\mathcal{C}_{2}}+J_{\mathcal{C}_{4}}-J_{\mathcal{C}_{5}}-J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}-J_{\mathcal{C}_{10}}+J_{\mathcal{C}_{11}}]\frac{\Delta\mu_{\rm
S}}{T_{0}}$ $\displaystyle=$ $\displaystyle J_{\rm
E}\Big{[}\frac{1}{T_{0}}-\frac{1}{(T_{0}+\delta T)}\Big{]}+J_{\rm
S}\Big{[}\frac{\mu_{\rm L\downarrow}}{T_{0}}-\frac{\mu_{\rm
L\uparrow}}{T_{0}}\Big{]}$ $\displaystyle=$ $\displaystyle J_{\rm
E}\mathcal{X}_{\rm E}+J_{\rm S}\mathcal{X}_{\rm S},$

where $\mathcal{X}_{\rm E}$ and $\mathcal{X}_{\rm S}$ are identified as
conjugate forces corresponding to the energy current $J_{\rm E}$, and spin
current $J_{\rm S}$, respectively. If we compare Eq. (47) and (63), we observe
that in both cases, the entropy production rate $\dot{\sigma}$ is the product
of the fluxes and forces: In Eq.(47), $\dot{\sigma}$ is in terms of the
microscopic fluxes ($J_{\mathcal{C}}$) and its conjugate forces
($\mathcal{X}_{\mathcal{C}}$), i.e., cycle affinities whereas in Eq.(63),
$\dot{\sigma}$ is in terms of the macroscopic fluxes (like the flow of spin
and energy current) and the associated phenomenological forces
($\mathcal{X}_{\rm E}$ and $\mathcal{X}_{\rm S}$). This is one of our central
results, showcasing the recovery of the phenomenological thermodynamic law of
entropy production in terms of generalized thermodynamic forces and fluxes
derived from the microscopic dynamical framework of the master equation
employing network cycle flux and forces. For small external bias $\delta T$
and $\Delta\mu_{\rm S}$, we can simplify Eq. (63) in the following form

$\displaystyle\dot{\sigma}\approx J_{\rm E}\left(\frac{\delta
T}{{T_{0}}^{2}}\right)+J_{\rm S}\left(\frac{\Delta\mu_{\rm S}}{T_{0}}\right),$
(64)

which is in accordance with the linear dependence of the entropy production
rate on the relevant thermodynamic forces. Similarly, the spin and the energy
currents within the linear response regime can be approximated as follows:

$\displaystyle J_{\rm S}$ $\displaystyle\approx$ $\displaystyle{\rm
U}(-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}})\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)}+(J^{-}_{\mathcal{C}_{1}}+J^{-}_{\mathcal{C}_{2}}+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}+J^{+}_{\mathcal{C}_{8}}+J^{-}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{11}})\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)},$ (65) $\displaystyle J_{\rm E}$
$\displaystyle\approx$ $\displaystyle{\rm
U^{2}}(J^{-}_{\mathcal{C}_{3}}+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{6}}+J^{-}_{\mathcal{C}_{7}}+J^{+}_{\mathcal{C}_{8}})\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)}+{\rm
U}(-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}})\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}.$ (66)

Comparing Eqs. (65) and (66) with the phenomenological linear law of
irreversible thermodynamics

$\displaystyle J^{\rm ph}_{\rm S}=L_{\rm SE}\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)}+L_{\rm SS}\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)},$ (67) $\displaystyle J^{\rm ph}_{\rm E}=L_{\rm
EE}\Bigg{(}\frac{\delta T}{k_{B}{T_{0}}^{2}}\Bigg{)}+L_{\rm
ES}\Bigg{(}\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Bigg{)},$ (68)

we identify the Onsager transport coefficients ($L$’s) in terms of the
microscopic circuit fluxes obtained from the network theory

$\displaystyle L_{\rm SE}$ $\displaystyle=$ $\displaystyle{\rm
U}(-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}})\equiv
L_{\rm ES},$ (69) $\displaystyle L_{\rm SS}$ $\displaystyle=$
$\displaystyle(J^{-}_{\mathcal{C}_{1}}+J^{-}_{\mathcal{C}_{2}}+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}$
(70) $\displaystyle+$ $\displaystyle
J^{+}_{\mathcal{C}_{8}}+J^{-}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{11}}),$
$\displaystyle L_{\rm EE}$ $\displaystyle=$ $\displaystyle{\rm
U^{2}}(J^{-}_{\mathcal{C}_{3}}+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{6}}+J^{-}_{\mathcal{C}_{7}}+J^{+}_{\mathcal{C}_{8}}).$
(71)

Equation (69) encapsulates the essence of the Onsager reciprocity relation.
Here we derive this relation by applying the quantum kinetic Pauli master
equation within the framework of network theory. It reveals that the BMS
quantum master equation is not a mere description of the dissipative dynamics
of the open quantum system; rather, it reproduces the reciprocity relation of
the linear law of irreversible thermodynamics, which obeys due to the time-
reversal symmetry of the stationary fluctuations. Here instead, it follows
from the properties of the network circuit fluxes between the forward
(counterclockwise) and reverse (clockwise) cycle trajectories.

Now, our aim is to establish the relationship between the coefficients of the
spin-Seebeck and the spin-Peltier effects. Under the zero spin current
condition, i.e., $J_{\rm S}=0$, we obtain from Eq. (67)

$\displaystyle L_{\rm SE}\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)}+L_{\rm SS}\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}=0,$ $\displaystyle{\rm or},\quad\kappa$
$\displaystyle\equiv\Bigg{(}\frac{\Delta\mu_{\rm S}}{\delta T}\Bigg{)}_{J_{\rm
S}=0}=-\frac{1}{T_{0}}\Bigg{(}\frac{L_{\rm SE}}{L_{\rm SS}}\Bigg{)}$ (72)
$\displaystyle\kappa$ $\displaystyle=-\frac{\rm
U}{T_{0}}\Big{(}-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}}\Big{)}/\Big{(}J^{-}_{\mathcal{C}_{1}}+J^{-}_{\mathcal{C}_{2}}$
(73)
$\displaystyle+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}+J^{+}_{\mathcal{C}_{8}}+J^{-}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{11}}\Big{)}.$

Here $\kappa=({\Delta\mu_{\rm S}}/{\delta T})_{J_{\rm S}=0}$ is the spin-
Seebeck coefficient or spin-thermoelectric power, defined as the change in the
spin bias voltage per unit change of temperature. Similarly, we may define the
spin-Peltier coefficient as

$\displaystyle\vartheta$ $\displaystyle=-\Bigg{(}\frac{J_{\rm E}}{J_{\rm
S}}\Bigg{)}_{{\delta T}=0}=-\frac{L_{\rm ES}}{L_{\rm SS}},$ (75)
$\displaystyle=-{\rm
U}\Big{(}-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}}\Big{)}/\Big{(}J^{-}_{\mathcal{C}_{1}}+J^{-}_{\mathcal{C}_{2}}$
$\displaystyle+J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{7}}+J^{+}_{\mathcal{C}_{8}}+J^{-}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{11}}\Big{)}.$

On the face of it, both Eqs. (72) and (75), appear exactly the same as their
classical counterparts, although their basis is completely different in
classical and quantum cases. Using Eqs. (73) and (75), we immediately conclude
that the classic Kelvin relation

$\displaystyle T_{0}\Bigg{(}\frac{\Delta\mu_{\rm S}}{\delta T}\Bigg{)}_{J_{\rm
S}=0}$ $\displaystyle=-\Bigg{(}\frac{J_{\rm E}}{J_{\rm S}}\Bigg{)}_{{\delta
T}=0}=-\frac{L_{\rm ES}}{L_{\rm SS}}=-\frac{L_{\rm SE}}{L_{\rm SS}},$ ${\rm
or},\quad T_{0}\kappa=\vartheta.$ (76)

equally holds for quantum thermocouples, connecting the two thermoelectric
effects, namely SSE and SPE. This is the hallmark of thermodynamics with its
universal generality. The generality prevails in the sense that all the
thermodynamic relations retain their forms in both classical and quantum
settings, with the only variation being in specific expressions that are used
to articulate them.

##  IV Operational Principles

To this end, we delve into the operational principles underlying the spin-
Seebeck and spin-Peltier effects.

###  IV.1 Spin-Seebeck effect

We observe the SSE when there is no spin bias voltage $\Delta\mu_{\rm S}=0$,
and a spin current is generated due to a temperature difference $\delta T$
between the upper and the lower terminals of the device. For $\Delta\mu_{\rm
S}=0$, Eq.(III.2) reduces to

$\displaystyle J_{\rm S}$ $\displaystyle=$
$\displaystyle-(J^{+}_{\mathcal{C}_{4}}-J^{-}_{\mathcal{C}_{4}})+(J^{+}_{\mathcal{C}_{5}}-J^{-}_{\mathcal{C}_{5}})+(J^{+}_{\mathcal{C}_{7}}-J^{-}_{\mathcal{C}_{7}})$
(77) $\displaystyle-$
$\displaystyle(J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}})$
$\displaystyle=$ $\displaystyle-
J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{7}}-J_{\mathcal{C}_{8}},$

while the $J_{\rm E}$ is still governed by Eq. (62). As a result, we identify
that cycle fluxes corresponding to subcycles $\mathcal{C}_{3}$,
$\mathcal{C}_{4}$, $\mathcal{C}_{5}$, $\mathcal{C}_{6}$, $\mathcal{C}_{7}$,
and $\mathcal{C}_{8}$ contribute to $J_{\rm E}$, while $\mathcal{C}_{4}$,
$\mathcal{C}_{5}$, $\mathcal{C}_{7}$, and $\mathcal{C}_{8}$ facilitate $J_{\rm
S}$. All the contributing cycle fluxes, energy, and spin currents, in
dimensionless units, along with all six microstate populations, are plotted in
Fig. 4 w.r.t the dimensionless temperature gradient $\delta T$. Although all
six cycles assist $J_{\rm E}$, the primary contribution comes from
$\mathcal{C}_{3}$, classified as the highest-rank cycle with a nonzero
contribution. In contrast, all the cycles appearing in $J_{\rm S}$ have equal
magnitudes but are lower in rank compared to $\mathcal{C}_{3}$ [Fig. 4a].
Consequently, the dimensionless spin current $J_{\rm S}$ observed in the SSE,
which is a linear combination of the four contributing cycles, is two orders
of magnitude less than the dimensionless energy current $J_{\rm E}$ [Figs. 4b
and 4c]. Upon setting $\Delta\mu_{\rm S}=0$ in Eqs. (65) and (66), approximate
expressions for $J_{\rm S}$ and $J_{\rm E}$ in the linear response regime are
reduced to

$J_{\rm S}\approx L_{\rm SE}\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)},$ (78)

and

$J_{\rm E}\approx L_{\rm EE}\Bigg{(}\frac{\delta
T}{k_{B}{T_{0}}^{2}}\Bigg{)},$ (79)

respectively, where $L_{\rm SE}$ and $L_{\rm EE}$ are given by Eqs. (69) and
(71) respectively. We observe that Eqs. (78) and (79) closely follow the
general Eqs. (77) and (62) [solid lines in Figs. 4b and 4c], but they start to
deviate [dash-dot lines] for large values of the temperature gradient.

Figure 4: Spin-Seebeck effect: All contributing (a) cycle fluxes (b)
populations (c) spin and (d) energy currents are plotted against dimensionless
thermal energy $k_{B}\delta T/\hbar\gamma$. The parameters used are as
follows: $\gamma_{\rm L}=\gamma_{\rm M}=\gamma_{\rm R}=\gamma$, ${\rm
U}=0.8\hbar\gamma$, $k_{B}T_{0}=4\hbar\gamma$,
$\varepsilon_{l\downarrow}=4.5\hbar\gamma$,
$\varepsilon_{l\uparrow}=0.5\hbar\gamma$, $\varepsilon_{u}=3\hbar\gamma$, and
$\mu_{\rm L\downarrow}=\mu_{\rm L\uparrow}=0$, $\mu_{\rm M}=0$.

Finally, we note that cycles ($\mathcal{C}_{4}$, $\mathcal{C}_{5}$,
$\mathcal{C}_{7}$, and $\mathcal{C}_{8}$) involving spin-flip processes
contribute to $J_{\rm S}$. For example, consider the dynamic steps of
$\mathcal{C}^{+}_{4}$: starting from the most populated state $\ket{00}$ as
shown in Fig.4(d), the system sequentially transitions to $\ket{0\uparrow}$
(where one spin-up electron tunnels from the left reservoir into the lower
QD), then to $\ket{1\uparrow}$ (where one electron tunnels from the middle
reservoir into the upper QD). The third step involves a spin-flip process
$\ket{1\uparrow}\rightarrow\ket{1\downarrow}$ by absorbing one magnon supplied
by the right reservoir. Subsequently, one spin-up electron tunnels into the
left reservoir ($\ket{1\uparrow}\rightarrow\ket{10}$), and finally, the system
returns to its initial state $\ket{00}$ by releasing one electron to the
middle reservoir. Thus, at the end of the full cycle, an integer spin-1 is
transferred from the left spinful electron reservoir to the right magnon bath.
Similarly, the clockwise circuit $\mathcal{C}^{-}_{4}$ represents the reverse
process, and both cycle trajectories additively contribute to the spin current
expression [Eq.(77)] in the SSE. The same holds true for other contributing
cycle trajectories mentioned in Eq. (77).

###  IV.2 Spin-Peltier effect

We observe the SPE in a scenario when $\delta T=0$ and an energy current is
generated due to a non-zero spin bias voltage. Putting $\delta T=0$ in Eq
(62), we obtain

$\displaystyle J_{\rm E}$ $\displaystyle=$ $\displaystyle{\rm
U}[(J^{+}_{\mathcal{C}_{4}}-J^{-}_{\mathcal{C}_{4}})+(J^{+}_{\mathcal{C}_{5}}-J^{-}_{\mathcal{C}_{5}})+(J^{+}_{\mathcal{C}_{7}}-J^{-}_{\mathcal{C}_{7}})$
(80) $\displaystyle+$
$\displaystyle(J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}})]$
$\displaystyle=$ $\displaystyle{\rm
U}[J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}],$
(81)

and the same Eq.(III.2) can be used to calculate $J_{\rm S}$.

Figure 5: Spin-Peltier effect: All contributing (a) cycle fluxes (b)
populations (c) spin and (d) energy currents are plotted against dimensionless
thermal energy $\Delta\mu_{\rm S}/\hbar\gamma$. The parameters used are as
follows: $\gamma_{\rm L}=\gamma_{\rm M}=\gamma_{\rm R}=\gamma$, ${\rm
U}=0.8\hbar\gamma$, $k_{B}T_{0}=4\hbar\gamma$,
$\varepsilon_{l\downarrow}=4.5\hbar\gamma$,
$\varepsilon_{l\uparrow}=0.5\hbar\gamma$, $\varepsilon_{u}=3\hbar\gamma$, and
$\mu_{\rm L\uparrow}=2.5\hbar\gamma$, $\mu_{\rm L\downarrow}=\mu_{\rm
L\uparrow}+\Delta\mu_{\rm S}$, $\delta T=0$,$\mu_{\rm M}=0$.

As a result, cycle fluxes corresponding to cycles $\mathcal{C}_{1}$,
$\mathcal{C}_{2}$, $\mathcal{C}_{4}$, $\mathcal{C}_{5}$, $\mathcal{C}_{7}$,
$\mathcal{C}_{8}$, $\mathcal{C}_{10}$, and $\mathcal{C}_{11}$ attribute to the
spin current, while only four cycles $\mathcal{C}_{4}$, $\mathcal{C}_{5}$,
$\mathcal{C}_{7}$, and $\mathcal{C}_{8}$ contribute to the energy current. All
supporting cycle fluxes, spin and energy currents in dimensionless units, and
the population of each eigenstate are plotted in Fig.5 against the
dimensionless spin bias voltage $\Delta\mu_{\rm S}$. In this case, the major
contribution to $J_{\rm S}$ is coming from cycle $\mathcal{C}_{1}$ which is
classified as the top-ranked cycle with a nonzero contribution. On the other
hand, all the cycles contributing to $J_{\rm E}$ are lower ranked cycles
relative to $\mathcal{C}_{1}$ with comparable magnitudes [Fig.5a].
Consequently, the dimensionless energy current $J_{\rm E}$ in SPE, is two
orders of magnitude less than the dimensionless spin current $J_{\rm S}$
[Fig.5b and 5c]. Upon substituting $\delta T=0$ in Eqs.(65) and (66),
approximate expressions for $J_{\rm E}$ and $J_{\rm S}$ in the linear response
regime take the form of

$J_{\rm E}\approx L_{\rm ES}\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)},$ (82)

and

$J_{\rm S}\approx L_{\rm SS}\Bigg{(}\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)},$ (83)

where $L_{\rm ES}$ and $L_{\rm SS}$ are given by Eqs. (69) and (70),
respectively. Similar to SSE, we note that for small bias, Eqs. (82) and (83)
closely follow the general Eqs. (81) and (61) [solid lines in Figs. 5b and
5c]. As mentioned earlier, each cycle $\mathcal{C}$ represents two paired
circuits, i.e., $\mathcal{C}^{+}$ and $\mathcal{C}^{-}$. In Fig.5a,
$J_{\mathcal{C}}<0$ implies that the flux corresponding to the
counterclockwise circuit $\mathcal{C}^{+}$ ($J^{+}_{\mathcal{C}}$) is less
than the flux corresponding to clockwise circuit $\mathcal{C}^{-}$
($J^{-}_{\mathcal{C}}$). Finally, we emphasize that it is the same set of four
cycles ($\mathcal{C}_{4}$, $\mathcal{C}_{5}$, $\mathcal{C}_{7}$, and
$\mathcal{C}_{8}$) that not only produces a weak spin current in SSE but also
accounts for generating a weak energy current in SPE.

###  IV.3 SSE and SPE: As thermodynamic cross-effect

In Section III.2, we have identified cycles
${\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{10},\mathcal{C}_{11}}$, which
can contribute solely to the spin current and not to energy current.
Conversely, cycles ${\mathcal{C}_{3},\mathcal{C}_{6}}$ are found to contribute
exclusively to the energy current and not to the spin current. Meanwhile,
cycles ${\mathcal{C}_{4},\mathcal{C}_{5},\mathcal{C}_{7},\mathcal{C}_{8}}$
have the quality to contribute to both energy as well as spin currents. To
gain a deeper understanding, it is essential to analyze the complete topology
of the network. Notably, we observe that all cycles contributing to $J_{\rm
S}$ must involve a spin-flip process, either
$\ket{0\uparrow}\leftrightarrow\ket{0\downarrow}$ or
$\ket{1\uparrow}\leftrightarrow\ket{1\downarrow}$, corresponding to the edges
$\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{4}}$ or
$\ket{\mathbb{5}}\leftrightarrow\ket{\mathbb{6}}$, respectively. Similarly,
cycles contributing to $J_{\rm E}$ must include the edges
$\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{5}}$
($\ket{0\uparrow}\leftrightarrow\ket{1\uparrow}$) or
$\ket{\mathbb{4}}\leftrightarrow\ket{\mathbb{6}}$
($\ket{0\downarrow}\leftrightarrow\ket{1\downarrow}$), enabling Coulomb
interaction between the upper and lower dots. This insight sheds light on why
the expressions for spin and energy currents, derived from the master
equation, take their particular forms [Cf. Eqs. (II) and (20)].

Figure 6: Schematic diagram of four subcycles
$\\{\mathcal{C}_{4},\mathcal{C}_{5},\mathcal{C}_{7},\mathcal{C}_{8}\\}$ which
are truly responsible for the spin-thermolectric cross-effect of SSE and SPE.
Special edges are marked in orange.

At this juncture, it’s imperative to underscore that cycles
${\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{9},\mathcal{C}_{10},\text{ and
}\mathcal{C}_{11}}$ share the edge
$\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{4}}$ or
$\ket{\mathbb{5}}\leftrightarrow\ket{\mathbb{6}}$, yet they do not contribute
to the spin current in the SSE due to zero cycle affinity. This results from
the fact that the circuit fluxes associated with the cycle trajectories
($\mathcal{C}^{+}$ and $\mathcal{C}^{-}$) are identical in both directions,
yielding a zero cycle flux. The same holds true for the SPE with cycles
$\mathcal{C}_{3}$, $\mathcal{C}_{6}$, and $\mathcal{C}_{9}$ in the absence of
a spin bias voltage, despite having the required edges. Intriguingly, cycle
$\mathcal{C}_{9}$ possesses both spin-flip and Coulomb-interaction edges [Fig.
2], yet it yields zero cycle fluxes due to its zero cycle affinity. Therefore,
the asymmetry in cycle affinity emerges as the primary thermodynamic driving
force and the foremost criterion for obtaining a nonzero cycle flux.
Conversely, cycles
${\mathcal{C}_{4},\mathcal{C}_{5},\mathcal{C}_{7},\mathcal{C}_{8}}$ exhibit
nonzero cycle affinity either in the absence of a temperature bias (SPE) or in
the absence of a spin bias voltage (SSE). Consequently, these four cycles
stand as the sole contributors to both SSE and SPE, featuring finite spin and
energy currents. This is attributed to their possession of both the spin-flip
edge ($\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{4}}$ or
$\ket{\mathbb{5}}\leftrightarrow\ket{\mathbb{6}}$) and the Coulomb-interaction
edge ($\ket{\mathbb{3}}\leftrightarrow\ket{\mathbb{5}}$ or
$\ket{\mathbb{4}}\leftrightarrow\ket{\mathbb{6}}$). To summarize, we affirm
that ${\mathcal{C}_{4},\mathcal{C}_{5},\mathcal{C}_{7},\mathcal{C}_{8}}$ stand
as the four pivotal cycles solely responsible for materializing the
thermodynamic cross-effect in the form of the spin-thermoelectric effect
within our simple minimal model of the quantum thermocouple. The consequences
of the interference effect between these major contributing cycles open
intriguing avenues, and future research directions could delve into the impact
of quantum coherence and entanglement effect Rao  _et al._ (2020); Whitney
(2016) on the device performance from the perspective of the network theory.

##  V Conclusions

The key findings of our present analysis are outlined as follows: (i) We
present a simple model of a quantum thermocouple that exhibits spin-
caloritronic effects based on three-terminal ultra-strong Coulomb-coupled
quantum dots. In contrast to four-terminal models, this minimal model mimics
both spin-dependent Seebeck and Peltier effects in complete analogy to
classical thermocouples, used to describe thermoelectric effects. In the
quantum case, distinct statistical properties of the thermal reservoirs play a
role akin to dissimilar metals in traditional thermocouples. (ii) We find out
that the expressions for spin and energy currents, derived from the Lindblad
master equation, completely agree with network theoretical results. However,
the quantum kinetic Pauli master equation serves as the basis for constructing
the thermodynamic network, encompassing joint system microstates and
associated transition rates. This is in stark contrast to classical network
theory, where microstates often result from coarse-graining procedures. Here
instead, they naturally emerge as eigenstates of the coupled quantum systems,
derived from the microscopic Hamiltonian description of the composite quantum
system. (iii) Benefiting from the generalized matrix tree theorem of the
algebraic graph, we not only unveil the fundamental operational principles
behind spin-Seebeck and spin-Peltier effects but also confirm the
applicability of well-known thermodynamic relations in nano-thermoelectric
devices. The validity of Onsager reciprocity and Kelvin relations for
thermoelectric coefficients underscore the universal generality of
thermodynamic principles in both classical and quantum realms. In the present
case, the above relations stem from the characteristic properties of forward
and backward cycle flux trajectories of the quantum thermodynamic network.
This is fundamentally different from the phenomenological classical laws of
irreversible thermodynamics that hinge on local equilibrium assumptions. (iv)
In this context, we stress the importance of network cycle flux, and cycle
affinity in establishing the macroscopic spin and energy currents in terms of
stochastic cycle currents. Cycle affinity, expressed as a ratio of transition
rates between forward and backward cycle trajectories, emerges as a
fundamental driving force behind nonzero cycle fluxes. Then, the cycle flux
ranking scheme powered by the microscopic or stochastic version of the entropy
production rate, sheds light on the origin of weak spin and energy currents in
spin-Seebeck and Peltier effects, respectively. (v) Finally, we identify four
non-intersecting cycles that are responsible for manifesting both reciprocal
effects of spin-thermoelectricity within our simple minimal model.
Characterized by special edges involving the spin-flip process and Coulomb
interaction between interacting quantum dots, these cycles pave the way for
underpinning the fundamental working principles of quantum thermocouples.

## Acknowledgements

We thank Sujan Kundu for the useful discussions. AG acknowledges financial
support from the Initiation grant of IITK (Grant No. IITK/CHM/2018513). N.G.
is thankful to CSIR for the fellowship. S.G. acknowledges the Ministry of
Education, Government of India, for the Prime Minister Research Fellowship
(PMRF).

## References

  * Rowe (1995) D. Rowe,  _CRC Handbook of Thermoelectrics_ (CRC-Press, 1995).
  * DiSalvo (1999) F. J. DiSalvo, Science 285, 703 (1999).
  * Goldsmid (2009) H. Goldsmid,  _Introduction to Thermoelectricity_, Springer Series in Materials Science (Springer Berlin Heidelberg, 2009).
  * Shakouri and Zebarjadi (2009) A. Shakouri and M. Zebarjadi, “Nanoengineered materials for thermoelectric energy conversion,” in  _Thermal Nanosystems and Nanomaterials_, edited by S. Volz (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009) pp. 225–299.
  * Dubi and Di Ventra (2011) Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 (2011).
  * Mazza  _et al._ (2015) F. Mazza, S. Valentini, R. Bosisio, G. Benenti, V. Giovannetti, R. Fazio, and F. Taddei, Phys. Rev. B 91, 245435 (2015).
  * Benenti  _et al._ (2017) G. Benenti, G. Casati, K. Saito, and R. Whitney, Physics Reports 694, 1 (2017), fundamental aspects of steady-state conversion of heat to work at the nanoscale.
  * Callen (1985) H. B. Callen,  _Thermodynamics and an Introduction to Thermostatistics_ , 2nd ed. (John Wiley & Sons, Inc., New York, 1985).
  * Landi and Paternostro (2021) G. T. Landi and M. Paternostro, Rev. Mod. Phys. 93, 035008 (2021).
  * Onsager (1931a) L. Onsager, Phys. Rev. 37, 405 (1931a).
  * Onsager (1931b) L. Onsager, Phys. Rev. 38, 2265 (1931b).
  * Callen (1948) H. B. Callen, Phys. Rev. 73, 1349 (1948).
  * Uchida  _et al._ (2008) K. Uchida, S. Takahashi, K. Harii, J. Ieda, W. Koshibae, K. Ando, S. Maekawa, and E. Saitoh, Nature 455, 778 (2008).
  * Wu  _et al._ (2015) S. M. Wu, J. E. Pearson, and A. Bhattacharya, Phys. Rev. Lett. 114, 186602 (2015).
  * Zhou  _et al._ (2021) W. Zhou, K. Yamamoto, A. Miura, R. Iguchi, Y. Miura, K.-i. Uchida, and Y. Sakuraba, Nature Materials 20, 463 (2021).
  * Flipse  _et al._ (2014) J. Flipse, F. K. Dejene, D. Wagenaar, G. E. W. Bauer, J. B. Youssef, and B. J. van Wees, Phys. Rev. Lett. 113, 027601 (2014).
  * Daimon  _et al._ (2016) S. Daimon, R. Iguchi, T. Hioki, E. Saitoh, and K.-i. Uchida, Nature Communications 7, 13754 (2016).
  * Ohnuma  _et al._ (2017) Y. Ohnuma, M. Matsuo, and S. Maekawa, Phys. Rev. B 96, 134412 (2017).
  * Bauer  _et al._ (2012) G. E. W. Bauer, E. Saitoh, and B. J. van Wees, Nature Materials 11, 391 (2012).
  * Boona  _et al._ (2014) S. R. Boona, R. C. Myers, and J. P. Heremans, Energy Environ. Sci. 7, 885 (2014).
  * Ronetti  _et al._ (2016) F. Ronetti, L. Vannucci, G. Dolcetto, M. Carrega, and M. Sassetti, Phys. Rev. B 93, 165414 (2016).
  * ichi UCHIDA (2021) K. ichi UCHIDA, Proceedings of the Japan Academy, Series B 97, 69 (2021).
  * Di Ventra (2008) M. Di Ventra,  _Electrical Transport in Nanoscale Systems_ (Cambridge University Press, 2008).
  * Nazarov and Blanter (2009) Y. Nazarov and Y. Blanter,  _Quantum Transport: Introduction to Nanoscience_ (Cambridge University Press, 2009).
  * Ihn (2010) T. Ihn,  _Semiconductor Nanostructures: Quantum States and Electronic Transport_ (OUP Oxford, 2010).
  * Heikkilä (2013) T. Heikkilä,  _The Physics of Nanoelectronics: Transport and Fluctuation Phenomena at Low Temperatures_, Oxford Master Series in Physics (OUP Oxford, 2013).
  * Ren (2013) J. Ren, Phys. Rev. B 88, 220406 (2013).
  * Whitney  _et al._ (2016) R. S. Whitney, R. Sánchez, F. Haupt, and J. Splettstoesser, Physica E: Low-dimensional Systems and Nanostructures 75, 257 (2016).
  * Whitney  _et al._ (2018) R. S. Whitney, R. Sánchez, and J. Splettstoesser, “Quantum thermodynamics of nanoscale thermoelectrics and electronic devices,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 175–206.
  * Wang  _et al._ (2022) L. Wang, Z. Wang, C. Wang, and J. Ren, Phys. Rev. Lett. 128, 067701 (2022).
  * van Houten  _et al._ (1992) H. van Houten, L. W. Molenkamp, C. W. J. Beenakker, and C. T. Foxon, Semiconductor Science and Technology 7, B215 (1992).
  * Lee  _et al._ (2016) M.-J. Lee, J.-H. Ahn, J. H. Sung, H. Heo, S. G. Jeon, W. Lee, J. Y. Song, K.-H. Hong, B. Choi, S.-H. Lee, and M.-H. Jo, Nature Communications 7, 12011 (2016).
  * Svilans  _et al._ (2016) A. Svilans, M. Leijnse, and H. Linke, Comptes Rendus Physique 17, 1096 (2016), mesoscopic thermoelectric phenomena / Phénomènes thermoélectriques mésoscopiques.
  * Erlingsson  _et al._ (2017) S. I. Erlingsson, A. Manolescu, G. A. Nemnes, J. H. Bardarson, and D. Sanchez, Phys. Rev. Lett. 119, 036804 (2017).
  * Patel  _et al._ (2020) A. Patel, D. Singh, Y. Sonvane, P. B. Thakor, and R. Ahuja, ACS Applied Materials & Interfaces 12, 46212 (2020).
  * Han  _et al._ (2020) W. Han, S. Maekawa, and X.-C. Xie, Nature Materials 19, 139 (2020).
  * Yang  _et al._ (2023) G. Yang, L. Sang, C. Zhang, N. Ye, A. Hamilton, M. S. Fuhrer, and X. Wang, Nature Reviews Physics 5, 466 (2023).
  * Esposito  _et al._ (2009) M. Esposito, K. Lindenberg, and C. V. den Broeck, Europhysics Letters 85, 60010 (2009).
  * Nakpathomkun  _et al._ (2010) N. Nakpathomkun, H. Q. Xu, and H. Linke, Phys. Rev. B 82, 235428 (2010).
  * Donsa  _et al._ (2014) S. Donsa, S. Andergassen, and K. Held, Phys. Rev. B 89, 125103 (2014).
  * Sothmann  _et al._ (2014) B. Sothmann, R. Sánchez, and A. N. Jordan, Nanotechnology 26, 032001 (2014).
  * Erdman  _et al._ (2017) P. A. Erdman, F. Mazza, R. Bosisio, G. Benenti, R. Fazio, and F. Taddei, Phys. Rev. B 95, 245432 (2017).
  * Esposito  _et al._ (2012) M. Esposito, N. Kumar, K. Lindenberg, and C. Van den Broeck, Phys. Rev. E 85, 031117 (2012).
  * Thierschmann  _et al._ (2015) H. Thierschmann, R. Sánchez, B. Sothmann, F. Arnold, C. Heyn, W. Hansen, H. Buhmann, and L. W. Molenkamp, Nature Nanotechnology 10, 854 (2015).
  * Jiang  _et al._ (2015) J.-H. Jiang, M. Kulkarni, D. Segal, and Y. Imry, Phys. Rev. B 92, 045309 (2015).
  * Zhang  _et al._ (2017) Y. Zhang, X. Zhang, Z. Ye, G. Lin, and J. Chen, Applied Physics Letters 110, 153501 (2017).
  * Ghosh  _et al._ (2022) S. Ghosh, N. Gupt, and A. Ghosh, Entropy 24 (2022).
  * Datta (2005) S. Datta,  _Quantum Transport: Atom to Transistor_, Ke xue qian yan cong shu (Cambridge University Press, 2005).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Gelbwaser-Klimovsky  _et al._ (2015) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, Adv. At. Mol. Opt. Phys. 64, 329 (2015).
  * Joulain  _et al._ (2016) K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Ghosh  _et al._ (2017) A. Ghosh, C. L. Latune, L. Davidovich, and G. Kurizki, Proc. Natl. Acad. Sci. U.S.A. 114, 12156 (2017).
  * Potts (2019) P. P. Potts, “Introduction to quantum thermodynamics (lecture notes),” (2019), arXiv:1906.07439 [quant-ph] .
  * Gupt  _et al._ (2021) N. Gupt, S. Bhattacharyya, and A. Ghosh, Phys. Rev. E 104, 054130 (2021).
  * Gupt  _et al._ (2022) N. Gupt, S. Bhattacharyya, B. Das, S. Datta, V. Mukherjee, and A. Ghosh, Phys. Rev. E 106, 024110 (2022).
  * Gupt  _et al._ (2023) N. Gupt, S. Ghosh, and A. Ghosh, Phys. Rev. E 108, 034305 (2023).
  * Schnakenberg (1976) J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976).
  * Hill and Chen (1975) T. L. Hill and Y. D. Chen, Proceedings of the National Academy of Sciences 72, 1291 (1975).
  * Kohler and Vollmerhaus (1980) H.-H. Kohler and E. Vollmerhaus, Journal of Mathematical Biology 9, 275 (1980).
  * Ren (2017) J. Ren, Frontiers of Physics 12, 120505 (2017).
  * Dutta  _et al._ (2020) A. Dutta, G. M. Schütz, and D. Chowdhury, Phys. Rev. E 101, 032402 (2020).
  * Vandaele  _et al._ (2017) K. Vandaele, S. J. Watzman, B. Flebus, A. Prakash, Y. Zheng, S. R. Boona, and J. P. Heremans, Materials Today Physics 1, 39 (2017).
  * Sothmann and Büttiker (2012) B. Sothmann and M. Büttiker, Europhysics Letters 99, 27001 (2012).
  * Strasberg (2022) P. Strasberg,  _Quantum Stochastic Thermodynamics: Foundations and Selected Applications_ (Oxford University Press, 2022).
  * Werlang  _et al._ (2014) T. Werlang, M. Marchiori, M. Cornelio, and D. Valente, Physical Review E 89, 062109 (2014).
  * Sinha  _et al._ (2011a) S. S. Sinha, A. Ghosh, and D. S. Ray, Phys. Rev. E 84, 041113 (2011a).
  * Sinha  _et al._ (2011b) S. S. Sinha, A. Ghosh, and D. S. Ray, Phys. Rev. E 84, 031118 (2011b).
  * Tutte (2001) W. T. Tutte,  _Graph Theory_ , Vol. 21 (Cambridge University Press, 2001).
  * Balakrishnan and Ranganathan (2012) R. Balakrishnan and K. Ranganathan,  _A Textbook of Graph Theory_ (Springer Science & Business Media, 2012).
  * King and Altman (1956) E. L. King and C. Altman, The Journal of Physical Chemistry 60, 1375 (1956).
  * Wu  _et al._ (2012) J. Wu, F. Liu, J. Ma, R. J. Silbey, and J. Cao, J. Chem. Phys. 137, 174111 (2012).
  * Einax  _et al._ (2011) M. Einax, M. Dierl, and A. Nitzan, The Journal of Physical Chemistry C 115, 21396 (2011).
  * Einax and Nitzan (2014) M. Einax and A. Nitzan, The Journal of Physical Chemistry C 118, 27226 (2014).
  * Kirchhoff (1847) G. Kirchhoff, Poggendorffs Ann. Phys. Chem. 148, 497 (1847).
  * Hill and Kedem (1966) T. L. Hill and O. Kedem, Journal of Theoretical Biology 10, 399 (1966).
  * Keizer (1972) J. Keizer, Journal of Statistical Physics 6, 67 (1972).
  * Ohga  _et al._ (2023) N. Ohga, S. Ito, and A. Kolchinsky, Phys. Rev. Lett. 131, 077101 (2023).
  * Rao  _et al._ (2020) D. D. B. Rao, A. Ghosh, D. Gelbwaser-Klimovsky, N. Bar-Gill, and G. Kurizki, New Journal of Physics 22, 083035 (2020).
  * Whitney (2016) R. S. Whitney, Entropy 18 (2016).

##  Appendix A Derivation of the Lindblad Quantum Master Equation

The total Hamiltonian of the overall three-terminal setup is given by

$H=H_{\rm S}+H_{\rm B}+H_{\rm I},$ (84)

where $H_{\rm S}$, $H_{\rm B}$, and $H_{\rm I}$ are the total Hamiltonian of
the system, bath, and system-bath interaction, respectively. The interaction
Hamiltonian $H_{\rm I}$ is defined as $H_{\rm I}=H_{\rm IL}+H_{\rm IM}+H_{\rm
IR}$, wherein $H_{\rm IL(R)}$ denotes the interaction of the lower quantum dot
(${\rm QD}_{l}$) with the left (right) bath, and $H_{\rm IM}$ represents the
interaction of the upper quantum dot (${\rm QD}_{u}$) with the middle bath.
The interaction Hamiltonian $H_{\rm I\alpha}$ ($\alpha={\rm L,M,R}$) for each
$\alpha$-th bath is given by Wang  _et al._ (2022); Gupt  _et al._ (2023)

$\displaystyle H_{\rm IL}$ $\displaystyle=$ $\displaystyle H_{\rm
IL\uparrow}+H_{\rm IL\downarrow},\;H_{\rm IL\uparrow}=\hbar\sum_{k}(t_{{\rm
L}k}b^{\dagger}_{{\rm L}\uparrow k}d_{l\uparrow}+t^{*}_{{\rm
L}k}d^{\dagger}_{l\uparrow}b_{{\rm L}\uparrow k}),\;H_{\rm
IL\downarrow}=\hbar\sum_{k}(t_{{\rm L}k}b^{\dagger}_{{\rm L}\downarrow
k}d_{l\downarrow}+t^{*}_{{\rm L}k}d^{\dagger}_{l\downarrow}b_{{\rm
L}\downarrow k}),$ $\displaystyle H_{\rm IM}$ $\displaystyle=$
$\displaystyle\hbar\sum_{k}(t_{{\rm M}k}b^{\dagger}_{{\rm
M}k}d_{u}+t^{*}_{{\rm M}k}d^{\dagger}_{u}b_{{\rm M}k}),\;H_{\rm
IR}=\hbar\sum_{q}(g_{{\rm R}q}a^{\dagger}_{{\rm
R}q}d^{\dagger}_{l\uparrow}d_{l\downarrow}+g^{*}_{{\rm
R}q}d^{\dagger}_{l\downarrow}d_{l\uparrow}a_{{\rm R}q}).$ (85)

To formulate the master equation, we begin with the derivation by considering
the von Neumann equation applied to the total density matrix $\rho_{\rm T}$ of
the combined system and reservoirs in the interaction picture, as given in
Breuer and Petruccione (2007)

$\frac{d\rho_{\rm T}}{dt}=-\frac{i}{\hbar}[H_{\rm I}(t),\rho_{\rm T}(t)].$
(86)

Integrating Eq.(86) and tracing out the bath degrees of freedom, the master
equation in terms of the reduced density matrix $\rho$ of the coupled quantum
dot system under the Born-Markov approximation can be written as Breuer and
Petruccione (2007)

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\Tr_{B}\int^{\infty}_{0}ds[H_{\rm
I}(t),[H_{\rm I}(t-s),\rho_{\rm T}(t)]].$ (87)

In Eq.(87), $\rho(t)=\Tr_{\rm B}\\{\rho_{\rm T}\\}\equiv\Tr_{\rm
B}\\{\rho(t)\otimes\rho_{\rm B}\\}$ where $\rho_{\rm B}=\rho_{\rm
L\uparrow}\otimes\rho_{\rm L\downarrow}\otimes\rho_{\rm M}\otimes\rho_{\rm
R}$, and $\Tr_{\rm B}\equiv\Tr_{{\rm L\uparrow},{\rm L\downarrow},{\rm M},{\rm
R}}$ stands for the trace over each bath degrees of freedom. As a result, we
can rewrite Eq. (87) as Gupt  _et al._ (2022); Ghosh  _et al._ (2022); Gupt
_et al._ (2023)

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\Tr_{{\rm L\uparrow},{\rm
L\downarrow},{\rm M},{\rm R}}\int^{\infty}_{0}ds[H_{\rm I}(t),[H_{\rm
I}(t-s),\rho(t)\otimes\rho_{\rm L\uparrow}\otimes\rho_{\rm
L\downarrow}\otimes\rho_{\rm M}\otimes\rho_{\rm R}]].$ (88)

Using the following relations Ghosh  _et al._ (2022); Gupt  _et al._ (2023)

$\displaystyle\Tr_{\rm L\uparrow}(b_{{\rm L}\uparrow k}(t)\rho_{\rm
L\uparrow})$ $\displaystyle=$ $\displaystyle 0=\Tr_{\rm
L\uparrow}(b^{\dagger}_{{\rm L}\uparrow k}(t)\rho_{\rm
L\uparrow}),\quad\quad\Tr_{\rm L\uparrow}(b_{{\rm L}\uparrow k}(t)\rho_{\rm
L\uparrow})=0=\Tr_{\rm L\uparrow}(b^{\dagger}_{{\rm L}\uparrow k}(t)\rho_{\rm
L\uparrow})$ (89) $\displaystyle\Tr_{\rm M}(b_{{\rm M}k}(t)\rho_{\rm M})$
$\displaystyle=$ $\displaystyle 0=\Tr_{\rm M}(b^{\dagger}_{{\rm
M}k}(t)\rho_{\rm M}),\;\quad\quad\Tr_{\rm R}(a_{{\rm R}q}(t)\rho_{\rm
R})=0=\Tr_{\rm R}(a^{\dagger}_{{\rm R}q}(t)\rho_{\rm R}),$ (90)

one can simplify Eq.(88) as Gupt  _et al._ (2022); Ghosh  _et al._ (2022);
Gupt  _et al._ (2023)

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\sum_{\beta}\Tr_{{\rm L\uparrow},{\rm
L\downarrow},{\rm M},{\rm R}}\int^{\infty}_{0}ds[H_{\rm I\beta}(t),[H_{\rm
I\beta}(t-s),\rho(t)\otimes\rho_{\rm L\uparrow}\otimes\rho_{\rm
L\downarrow}\otimes\rho_{\rm M}\otimes\rho_{\rm R}]],\quad\beta={\rm
L\uparrow,L\downarrow,M,R}.$ (91)

Now, we use system operators in the interaction picture as

$\displaystyle d_{\rm i}(t)$ $\displaystyle=$ $\displaystyle e^{{iH_{\rm
S}t}/{\hbar}}d_{\rm i}e^{{-iH_{\rm
S}t}/{\hbar}}=\sum_{\\{\varepsilon_{\mathbb{ji}}\\}}e^{{-i{\varepsilon_{\mathbb{ji}}}t}/{\hbar}}d_{\rm
i},$ $\displaystyle d^{\dagger}_{\rm i}(t)$ $\displaystyle=$ $\displaystyle
e^{{iH_{\rm S}t}/{\hbar}}d^{\dagger}_{\rm i}e^{{-iH_{\rm
S}t}/{\hbar}}=\sum_{\\{\varepsilon_{\mathbb{ji}}\\}}e^{{i{\varepsilon_{\mathbb{ji}}}t}/{\hbar}}d^{\dagger}_{\rm
i},\quad{\rm i}=l\uparrow,l\downarrow,u$ (92)

where
$\varepsilon_{\mathbb{ji}}=\varepsilon_{\mathbb{j}}-\varepsilon_{\mathbb{i}}>0$
is the energy required for the transition between state $|\mathbb{i}\rangle$
and $|\mathbb{j}\rangle$ driven by their respective bath. Similarly, one can
write the expressions for the bath operators in the interaction picture. With
all these given prescriptions, we have simplified the Eq.(91), resulting in
the Lindblad form of the quantum master equation as follows:

$\frac{d\rho}{dt}=\mathcal{L}_{\rm L\uparrow}[\rho]+\mathcal{L}_{\rm
L\downarrow}[\rho]+\mathcal{L}_{\rm M}[\rho]+\mathcal{L}_{\rm R}[\rho],$ (93)

The explicit forms of the Lindblad super operator $\mathcal{L}$ in the above
equation are given by

$\displaystyle\mathcal{L}_{\rm L\uparrow}[\rho]=\sum_{\\{\varepsilon_{\rm
L\uparrow}\\}}\gamma_{\rm L}\Big{[}f(\varepsilon_{\rm L\uparrow},\mu_{\rm
L\uparrow},T_{\rm L})\Big{(}d^{\dagger}_{l\uparrow}(\varepsilon_{\rm
L\uparrow})\rho d_{l\uparrow}(\varepsilon_{\rm
L\uparrow})-\frac{1}{2}\\{d_{l\uparrow}(\varepsilon_{\rm
L\uparrow})d^{\dagger}_{l\uparrow}(\varepsilon_{\rm
L\uparrow}),\rho\\}\Big{)}$ $\displaystyle+(1-f(\varepsilon_{\rm
L\uparrow},\mu_{\rm L\uparrow},T_{\rm
L}))\Big{(}d_{l\uparrow}(\varepsilon_{\rm L\uparrow})\rho
d^{\dagger}_{l\uparrow}(\varepsilon_{\rm
L\uparrow})-\frac{1}{2}\\{d^{\dagger}_{l\uparrow}(\varepsilon_{\rm
L\uparrow})d_{l\uparrow}(\varepsilon_{\rm L\uparrow}),\rho\\}\Big{)}\Big{]},$
(94) $\displaystyle\mathcal{L}_{\rm
L\downarrow}[\rho]=\sum_{\\{\varepsilon_{\rm L\downarrow}\\}}\gamma_{\rm
L}\Big{[}f(\varepsilon_{\rm L\downarrow},\mu_{\rm L\downarrow},T_{\rm
L})\Big{(}d^{\dagger}_{l\downarrow}(\varepsilon_{\rm L\downarrow})\rho
d_{l\downarrow}(\varepsilon_{\rm
L\downarrow})-\frac{1}{2}\\{d_{l\downarrow}(\varepsilon_{\rm
L\downarrow})d^{\dagger}_{l\downarrow}(\varepsilon_{\rm
L\downarrow}),\rho\\}\Big{)}$ $\displaystyle+(1-f(\varepsilon_{\rm
L\downarrow},\mu_{\rm L\downarrow},T_{\rm
L}))\Big{(}d_{l\downarrow}(\varepsilon_{\rm L\downarrow})\rho
d^{\dagger}_{l\downarrow}(\varepsilon_{\rm
L\downarrow})-\frac{1}{2}\\{d^{\dagger}_{l\downarrow}(\varepsilon_{\rm
L\downarrow})d_{l\downarrow}(\varepsilon_{\rm
L\downarrow}),\rho\\}\Big{)}\Big{]},$ (95) $\displaystyle\mathcal{L}_{\rm
M}[\rho]=\sum_{\\{\varepsilon_{\rm M}\\}}\gamma_{\rm
M}\Big{[}f(\varepsilon_{\rm M},\mu_{\rm M},T_{\rm
L})\Big{(}d^{\dagger}_{u}(\varepsilon_{\rm M})\rho d_{u}(\varepsilon_{\rm
M})-\frac{1}{2}\\{d_{u}(\varepsilon_{\rm M})d^{\dagger}_{u}(\varepsilon_{\rm
M}),\rho\\}\Big{)}$ $\displaystyle+(1-f(\varepsilon_{\rm M},\mu_{\rm M},T_{\rm
M}))\Big{(}d_{u}(\varepsilon_{\rm M})\rho d^{\dagger}_{u}(\varepsilon_{\rm
M})-\frac{1}{2}\\{d^{\dagger}_{u}(\varepsilon_{\rm M})d_{u}(\varepsilon_{\rm
M}),\rho\\}\Big{)}\Big{]},$ (96) $\displaystyle\mathcal{L}_{\rm
R}[\rho]=\sum_{\\{\varepsilon_{\rm R}\\}}\gamma_{\rm
R}\Big{[}n(\varepsilon_{\rm R},T_{\rm
R})\Big{(}V^{\dagger}_{l}(\varepsilon_{\rm R})\rho V_{l}(\varepsilon_{\rm
R})-\frac{1}{2}\\{V_{l}(\varepsilon_{\rm R})V^{\dagger}_{l}(\varepsilon_{\rm
R}),\rho\\}\Big{)}$ $\displaystyle+(n(\varepsilon_{\rm R},T_{\rm
R})+1)\Big{(}V_{l}(\varepsilon_{\rm R})\rho V^{\dagger}_{l}(\varepsilon_{\rm
R})-\frac{1}{2}\\{V^{\dagger}_{l}(\varepsilon_{\rm R})V_{l}(\varepsilon_{\rm
R}),\rho\\}\Big{)}\Big{]},$ (97)

where the operators $V_{l}=d^{\dagger}_{l\uparrow}d_{l\downarrow}$ and
$V^{\dagger}_{l}=d^{\dagger}_{l\downarrow}d_{l\uparrow}$ are responsible for
the transition between spin-up ($\uparrow$) and spin-down ($\downarrow$)
states. The transition rates corresponding to their respective bath are
characterized by the various $\gamma$’s. The explicit form of all $\gamma$’s
in terms of system-bath coupling constants can be calculated by Fermi’s golden
rule, as $\gamma_{\rm L}=2\pi\hbar\sum_{k}|t_{{\rm
L}k}|^{2}\delta\big{(}\varepsilon-\epsilon_{{\rm L\sigma}k}\big{)}$, where
$\sigma=\\{\uparrow,\downarrow\\}$, $\gamma_{\rm M}=2\pi\hbar\sum_{k}|t_{{\rm
M}k}|^{2}\delta\big{(}\varepsilon-\epsilon_{{\rm M}k}\big{)}$, and
$\gamma_{\rm R}=2\pi\hbar\sum_{q}|g_{{\rm
R}q}|^{2}\delta\big{(}\varepsilon-\epsilon_{{\rm R}q}\big{)}$ Gupt  _et al._
(2023). The functions
$f(\varepsilon,\mu,T)=[e^{(\varepsilon-\mu)/k_{B}T}+1]^{-1}$ is the Fermi-
Dirac distribution function corresponding to the left (L) and middle (M) bath
with the transition energy $\varepsilon$, chemical potential $\mu$, and
equilibrium bath temperature $T$. Similarly, the function
$n(\varepsilon,T)=[e^{\varepsilon/k_{B}T}-1]^{-1}$ is the Bose-Einstein
distribution function corresponding to the right (R) bath with the transition
energy $\varepsilon$ and reservoir temperature $T$. The distribution functions
are defined as the bath correlation functions and can be calculated as
$\langle{b^{\dagger}b}\rangle=\Tr_{\rm L\sigma(M)}({b^{\dagger}b\rho_{\rm
L\sigma(M)}})=f_{\rm L\sigma(M)}$ and $\langle{bb^{\dagger}}\rangle=\Tr_{\rm
L\sigma(M)}({bb^{\dagger}\rho_{\rm L\sigma(M)}})=1-f_{\rm L\sigma(M)}$ for the
left (middle) bath and $\langle{a^{\dagger}a}\rangle=\Tr_{\rm
R}({a^{\dagger}a\rho_{\rm R}})=n_{\rm R}$,
$\langle{aa^{\dagger}}\rangle=\Tr_{\rm R}({aa^{\dagger}\rho_{\rm R}})=1+n_{\rm
R}$ for the right bath Gupt  _et al._ (2023). The operators $b$ and
$b^{\dagger}$ follow anti-commutation relation whereas the operators $a$ and
$a^{\dagger}$ follow commutation relation, and $k_{B}$ is the Boltzmann
constant. The energies needed for the transitions which are driven by the left
and middle baths are $\varepsilon_{\rm
L\uparrow}=\\{\varepsilon_{\mathbb{31}},\varepsilon_{\mathbb{52}}\\}$,
$\varepsilon_{\rm
L\downarrow}=\\{\varepsilon_{\mathbb{41}},\varepsilon_{\mathbb{62}}\\}$, and
$\varepsilon_{\rm
M}=\\{\varepsilon_{\mathbb{21}},\varepsilon_{\mathbb{53}},\varepsilon_{\mathbb{64}}\\}$
respectively, while the energies required for the transitions triggered by the
right bath are $\varepsilon_{\rm
R}=\\{\varepsilon_{\mathbb{43}},\varepsilon_{\mathbb{65}}\\}$. Note that one
can express the various system creation and annihilation operators and their
combinations in the following forms $|\mathbb{i}\rangle\langle\mathbb{j}|$
($\mathbb{i}\neq\mathbb{j}$, $\mathbb{i},\;\mathbb{j}=\mathbb{1,2,3,4,5,6}$),
which are given by:

$\displaystyle d^{\dagger}_{l\uparrow}$ $\displaystyle=$
$\displaystyle|\mathbb{3}\rangle\langle\mathbb{1}|+|\mathbb{5}\rangle\langle\mathbb{2}|,\quad
d_{l\uparrow}=|\mathbb{1}\rangle\langle\mathbb{3}|+|\mathbb{2}\rangle\langle\mathbb{5}|,$
$\displaystyle d^{\dagger}_{l\downarrow}$ $\displaystyle=$
$\displaystyle|\mathbb{4}\rangle\langle\mathbb{1}|+|\mathbb{6}\rangle\langle\mathbb{2}|,\quad
d_{l\downarrow}=|\mathbb{1}\rangle\langle\mathbb{4}|+|\mathbb{2}\rangle\langle\mathbb{6}|,$
$\displaystyle d^{\dagger}_{u}$ $\displaystyle=$
$\displaystyle|\mathbb{2}\rangle\langle\mathbb{1}|+|\mathbb{5}\rangle\langle\mathbb{3}|+|\mathbb{6}\rangle\langle\mathbb{4}|,\quad
d_{u}=|\mathbb{1}\rangle\langle\mathbb{2}|+|\mathbb{3}\rangle\langle\mathbb{5}|+|\mathbb{4}\rangle\langle\mathbb{6}|,$
$\displaystyle V^{\dagger}_{l}$ $\displaystyle=$ $\displaystyle
d^{\dagger}_{l\downarrow}d_{l\uparrow}=|\mathbb{4}\rangle\langle\mathbb{3}|+|\mathbb{6}\rangle\langle\mathbb{5}|,\quad
V_{l}=d^{\dagger}_{l\uparrow}d_{l\downarrow}=|\mathbb{3}\rangle\langle\mathbb{4}|+|\mathbb{5}\rangle\langle\mathbb{6}|.$
(98)

Finally, Eqs. (9) to (14) in the main text can be derived with the help of
Eqs. (93)-(A) in the following manner. For example

$\frac{dP_{\mathbb{1}}}{dt}=\langle\mathbb{1}|\frac{d\rho}{dt}|\mathbb{1}\rangle=\langle\mathbb{1}|\mathcal{L}_{\rm
L\uparrow}[\rho]|\mathbb{1}\rangle+\langle\mathbb{1}|\mathcal{L}_{\rm
L\downarrow}[\rho]|\mathbb{1}\rangle+\langle\mathbb{1}|\mathcal{L}_{\rm
M}[\rho]|\mathbb{1}\rangle+\langle\mathbb{1}|\mathcal{L}_{\rm
R}[\rho]|\mathbb{1}\rangle,$ (99)

where the terms

$\displaystyle\langle\mathbb{1}|\mathcal{L}_{\rm
L\uparrow}[\rho]|\mathbb{1}\rangle$ $\displaystyle=$
$\displaystyle{\gamma_{\rm L}f(\varepsilon_{\mathbb{31}},\mu_{\rm
L\uparrow},T_{\rm
L})}\Big{(}-\frac{1}{2}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{1}|\rho|\mathbb{1}\rangle-\frac{1}{2}\langle\mathbb{1}|\rho|\mathbb{1}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}+\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{31}},\mu_{\rm L\uparrow},T_{\rm
L}))\Big{(}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{3}|\rho|\mathbb{3}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}$
(100) $\displaystyle=$ $\displaystyle\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{31}},\mu_{\rm L\uparrow},T_{\rm
L}))P_{\mathbb{3}}-\gamma_{L}f(\varepsilon_{\mathbb{31}},\mu_{\rm
L\uparrow},T_{\rm
L})P_{\mathbb{1}}=k_{\mathbb{13}}P_{\mathbb{3}}-k_{\mathbb{31}}P_{\mathbb{1}}\equiv
J_{\mathbb{13}},$ $\displaystyle\langle\mathbb{1}|\mathcal{L}_{\rm
L\downarrow}[\rho]|\mathbb{1}\rangle$ $\displaystyle=$
$\displaystyle{\gamma_{\rm L}f(\varepsilon_{\mathbb{41}},\mu_{\rm
L\downarrow},T_{\rm
L})}\Big{(}-\frac{1}{2}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{1}|\rho|\mathbb{1}\rangle-\frac{1}{2}\langle\mathbb{1}|\rho|\mathbb{1}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}+\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{41}},\mu_{\rm L\downarrow},T_{\rm
L}))\Big{(}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{4}|\rho|\mathbb{4}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}$
(101) $\displaystyle=$ $\displaystyle\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{41}},\mu_{\rm L\downarrow},T_{\rm
L}))P_{\mathbb{4}}-\gamma_{L}f(\varepsilon_{\mathbb{41}},\mu_{\rm
L\downarrow},T_{\rm
L})P_{\mathbb{1}}=k_{\mathbb{14}}P_{\mathbb{4}}-k_{\mathbb{41}}P_{\mathbb{1}}\equiv
J_{\mathbb{14}},$ $\displaystyle\langle\mathbb{1}|\mathcal{L}_{\rm
M}[\rho]|\mathbb{1}\rangle$ $\displaystyle=$ $\displaystyle{\gamma_{\rm
L}f(\varepsilon_{\mathbb{21}},\mu_{\rm L\downarrow},T_{\rm
L})}\Big{(}-\frac{1}{2}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{1}|\rho|\mathbb{1}\rangle-\frac{1}{2}\langle\mathbb{1}|\rho|\mathbb{1}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}+\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{21}},\mu_{\rm L\downarrow},T_{\rm
L}))\Big{(}\langle\mathbb{1}|\mathbb{1}\rangle\langle\mathbb{2}|\rho|\mathbb{2}\rangle\langle\mathbb{1}|\mathbb{1}\rangle\Big{)}$
(102) $\displaystyle=$ $\displaystyle\gamma_{\rm
L}(1-f(\varepsilon_{\mathbb{21}},\mu_{\rm L\downarrow},T_{\rm
L}))P_{\mathbb{2}}-\gamma_{L}f(\varepsilon_{\mathbb{21}},\mu_{\rm
L\downarrow},T_{\rm
L})P_{\mathbb{1}}=k_{\mathbb{12}}P_{\mathbb{2}}-k_{\mathbb{21}}P_{\mathbb{1}}\equiv
J_{\mathbb{12}},$ $\displaystyle\langle\mathbb{1}|\mathcal{L}_{\rm
R}[\rho]|\mathbb{1}\rangle$ $\displaystyle=$ $\displaystyle 0.$ (103)

Similarly, one can derive time evolution equations for the population of the
other $\mathbb{i}$-th states. Under the steady state, $dP_{\mathbb{i}}/dt=0$
($\mathbb{i=1,2,..,6}$), and we have

$\displaystyle\frac{dP_{\mathbb{1}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{13}}P_{\mathbb{3}}-k_{\mathbb{31}}P_{\mathbb{1}})+(k_{\mathbb{14}}P_{\mathbb{4}}-k_{\mathbb{41}}P_{\mathbb{1}})+(k_{\mathbb{12}}P_{\mathbb{2}}-k_{\mathbb{21}}P_{\mathbb{1}})=J_{\mathbb{13}}+J_{\mathbb{14}}+J_{\mathbb{12}}=0,$
(104) $\displaystyle\frac{dP_{\mathbb{2}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{25}}P_{\mathbb{5}}-k_{\mathbb{52}}P_{\mathbb{2}})+(k_{\mathbb{26}}P_{\mathbb{6}}-k_{\mathbb{62}}P_{\mathbb{2}})+(k_{\mathbb{21}}P_{\mathbb{1}}-k_{\mathbb{12}}P_{\mathbb{2}})=J_{\mathbb{25}}+J_{\mathbb{26}}+J_{\mathbb{21}}=0,$
(105) $\displaystyle\frac{dP_{\mathbb{3}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{31}}P_{\mathbb{1}}-k_{\mathbb{13}}P_{\mathbb{3}})+(k_{\mathbb{35}}P_{\mathbb{5}}-k_{\mathbb{53}}P_{\mathbb{3}})+(k_{\mathbb{34}}P_{\mathbb{4}}-k_{\mathbb{43}}P_{\mathbb{3}})=J_{\mathbb{31}}+J_{\mathbb{35}}+J_{\mathbb{34}}=0,$
(106) $\displaystyle\frac{dP_{\mathbb{4}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{41}}P_{\mathbb{1}}-k_{\mathbb{14}}P_{\mathbb{4}})+(k_{\mathbb{46}}P_{\mathbb{6}}-k_{\mathbb{64}}P_{\mathbb{4}})+(k_{\mathbb{43}}P_{\mathbb{3}}-k_{\mathbb{34}}P_{\mathbb{4}})=J_{\mathbb{41}}+J_{\mathbb{46}}+J_{\mathbb{43}}=0,$
(107) $\displaystyle\frac{dP_{\mathbb{5}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{52}}P_{\mathbb{2}}-k_{\mathbb{25}}P_{\mathbb{5}})+(k_{\mathbb{53}}P_{\mathbb{3}}-k_{\mathbb{35}}P_{\mathbb{5}})+(k_{\mathbb{56}}P_{\mathbb{6}}-k_{\mathbb{65}}P_{\mathbb{5}})=J_{\mathbb{52}}+J_{\mathbb{53}}+J_{\mathbb{56}}=0,$
(108) $\displaystyle\frac{dP_{\mathbb{6}}}{dt}$ $\displaystyle=$
$\displaystyle(k_{\mathbb{62}}P_{\mathbb{2}}-k_{\mathbb{26}}P_{\mathbb{6}})+(k_{\mathbb{64}}P_{\mathbb{4}}-k_{\mathbb{46}}P_{\mathbb{6}})+(k_{\mathbb{65}}P_{\mathbb{5}}-k_{\mathbb{56}}P_{\mathbb{6}})=J_{\mathbb{62}}+J_{\mathbb{64}}+J_{\mathbb{65}}=0.$
(109)

In the present case, there is no particle exchange between the quantum dots,
implying no particle current due to the middle reservoir. So, the steady-state
energy (heat) current through the middle reservoir within the Born-Markov-
Secular (BMS) master equation can be defined as Ghosh  _et al._ (2022); Wang
_et al._ (2022)

$J_{\rm E}=J^{\rm M}_{\rm E}=\Tr\\{\mathcal{L}_{\rm M}[\rho]H_{\rm S}\\},$
(110)

where the system Hamiltonian $H_{\rm S}$ has the following form, $H_{\rm
S}=\sum_{\mathbb{i}}\varepsilon_{\mathbb{i}}|\mathbb{i}\rangle\langle\mathbb{i}|$,
with $\varepsilon_{\mathbb{i}}$ being the energy of the $\mathbb{i}$-th state.
Using Eq.(110), the expression for t$J_{\rm E}$ can be calculated as

$\displaystyle J_{\rm
E}=\varepsilon_{\mathbb{21}}J_{\mathbb{21}}+\varepsilon_{\mathbb{53}}J_{\mathbb{53}}+\varepsilon_{\mathbb{64}}J_{\mathbb{64}}.$
(111)

At the steady state, one may verify
$J_{\mathbb{12}}=-J_{\mathbb{21}}=J_{\mathbb{53}}+J_{\mathbb{64}}$. As a
result, Eq. (111) reduce to Eq.(20) of the main text:

$\displaystyle J_{\rm E}$ $\displaystyle=$
$\displaystyle\varepsilon_{u}(-J_{\mathbb{12}})+(\varepsilon_{u}+{\rm
U})J_{\mathbb{53}}+(\varepsilon_{u}+{\rm U})J_{\mathbb{64}},$ (112)
$\displaystyle=$
$\displaystyle-\varepsilon_{u}J_{\mathbb{12}}+(\varepsilon_{u}+{\rm
U})(J_{\mathbb{53}}+J_{\mathbb{64}}),$ $\displaystyle=$ $\displaystyle{\rm
U}(J_{\mathbb{53}}+J_{\mathbb{64}}).$

Similarly, the steady-state spin current due to the left and right reservoirs
can be defined as Wang  _et al._ (2022):

$\displaystyle J^{\rm L}_{\rm S}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}\Tr\\{d^{\dagger}_{l\downarrow}d_{l\downarrow}\mathcal{L}_{\rm
L\downarrow}[\rho]\\}-\Tr\\{d^{\dagger}_{l\uparrow}d_{l\uparrow}\mathcal{L}_{\rm
L\uparrow}[\rho]\\}\Big{)},$ (113) $\displaystyle J^{\rm R}_{\rm S}$
$\displaystyle=$ $\displaystyle\Tr\\{V^{\dagger}_{l}V_{l}\mathcal{L}_{\rm
R}[\rho]\\}.$ (114)

Now, using Eq.(A), one can derive the expressions for $J^{\rm L}_{\rm S}$ and
$J^{\rm R}_{\rm S}$ from Eqs. (113) and (114) in the follwing forms:

$\displaystyle J^{\rm L}_{\rm S}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{[}(J_{\mathbb{41}}+J_{\mathbb{62}})-(J_{\mathbb{31}}+J_{\mathbb{52}})\Big{]},$
(115) $\displaystyle J^{\rm R}_{\rm S}$ $\displaystyle=$ $\displaystyle
J_{\mathbb{43}}+J_{\mathbb{65}}.$ (116)

In the steady-state,
$J_{\mathbb{41}}+J_{\mathbb{62}}=-(J_{\mathbb{31}}+J_{\mathbb{52}})=J_{\mathbb{34}}+J_{\mathbb{56}}$.
As a result, we get Eq.(II) of the main text as the steady-state spin current

$\displaystyle J_{\rm S}$ $\displaystyle=$ $\displaystyle J^{\rm L}_{\rm
S}=\frac{1}{2}\Big{[}(J_{\mathbb{41}}+J_{\mathbb{62}})+(J_{\mathbb{41}}+J_{\mathbb{62}})\Big{]}$
(117) $\displaystyle=$ $\displaystyle
J_{\mathbb{41}}+J_{\mathbb{62}}=J_{\mathbb{34}}+J_{\mathbb{56}}$
$\displaystyle=$ $\displaystyle-(J_{\mathbb{43}}+J_{\mathbb{65}})=-J^{\rm
R}_{\rm S}.$ (118)

##  Appendix B Derivation of the entropy production rate

From Eq.(46), the steady-state entropy production rate can be written as
follows Schnakenberg (1976); Landi and Paternostro (2021):

$\displaystyle\dot{\sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{2}k_{B}\sum_{\mathbb{i,j}}J_{\mathbb{ij}}\ln\Big(\frac{k_{\mathbb{ij}}}{k_{\mathbb{ij}}}\Big{missing})$
(119) $\displaystyle=$
$\displaystyle\frac{1}{2}k_{B}\Bigg{[}J_{\mathbb{31}}\ln\Big(\frac{k_{\mathbb{31}}}{k_{\mathbb{13}}}\Big{missing})+J_{\mathbb{13}}\ln\Big(\frac{k_{\mathbb{13}}}{k_{\mathbb{31}}}\Big{missing})+J_{\mathbb{41}}\ln\Big(\frac{k_{\mathbb{41}}}{k_{\mathbb{14}}}\Big{missing})+J_{\mathbb{14}}\ln\Big(\frac{k_{\mathbb{14}}}{k_{\mathbb{41}}}\Big{missing})+J_{\mathbb{43}}\ln\Big(\frac{k_{\mathbb{43}}}{k_{\mathbb{34}}}\Big{missing})+J_{\mathbb{34}}\ln\Big(\frac{k_{\mathbb{34}}}{k_{\mathbb{43}}}\Big{missing})+.....\Bigg{]}$

As we have mentioned in the main text $J_{\mathbb{ij}}=-J_{\mathbb{ji}}$ for
all $\mathbb{i}$ and $\mathbb{j}$ ($\mathbb{i\neq j}$), so $\dot{\sigma}$ will
be equal to

$\displaystyle\dot{\sigma}$ $\displaystyle=$ $\displaystyle
k_{B}\Bigg{[}J_{\mathbb{31}}\ln\Big(\frac{k_{\mathbb{31}}}{k_{\mathbb{13}}}\Big{missing})+J_{\mathbb{41}}\ln\Big(\frac{k_{\mathbb{41}}}{k_{\mathbb{14}}}\Big{missing})+J_{\mathbb{43}}\ln\Big(\frac{k_{\mathbb{43}}}{k_{\mathbb{34}}}\Big{missing})+J_{\mathbb{21}}\ln\Big(\frac{k_{\mathbb{21}}}{k_{\mathbb{12}}}\Big{missing})+J_{\mathbb{52}}\ln\Big(\frac{k_{\mathbb{52}}}{k_{\mathbb{25}}}\Big{missing})+J_{\mathbb{64}}\ln\Big(\frac{k_{\mathbb{64}}}{k_{\mathbb{46}}}\Big{missing})+.....\Bigg{]}$
(120) $\displaystyle=$ $\displaystyle
k_{B}\Bigg{[}(J^{+}_{\mathcal{C}_{1}}-J^{-}_{\mathcal{C}_{1}}+J^{+}_{\mathcal{C}_{3}}-J^{-}_{\mathcal{C}_{3}}+J^{+}_{\mathcal{C}_{4}}-J^{-}_{\mathcal{C}_{4}}+J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}}+J^{+}_{\mathcal{C}_{11}}-J^{-}_{\mathcal{C}_{11}})\ln\Big(\frac{k_{\mathbb{31}}}{k_{\mathbb{13}}}\Big{missing})+(-J^{+}_{\mathcal{C}_{1}}+J^{-}_{\mathcal{C}_{1}}+J^{+}_{\mathcal{C}_{5}}-J^{-}_{\mathcal{C}_{5}}+J^{+}_{\mathcal{C}_{6}}-J^{-}_{\mathcal{C}_{6}}$
$\displaystyle+$ $\displaystyle
J^{+}_{\mathcal{C}_{7}}-J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{11}}+J^{-}_{\mathcal{C}_{11}})\ln\Big(\frac{k_{\mathbb{41}}}{k_{\mathbb{14}}}\Big{missing})+(J^{+}_{\mathcal{C}_{1}}-J^{-}_{\mathcal{C}_{1}}-J^{+}_{\mathcal{C}_{5}}+J^{-}_{\mathcal{C}_{5}}+J^{+}_{\mathcal{C}_{8}}-J^{-}_{\mathcal{C}_{8}}-J^{+}_{\mathcal{C}_{9}}-J^{-}_{\mathcal{C}_{9}}-J^{+}_{\mathcal{C}_{10}}+J^{-}_{\mathcal{C}_{10}})\ln\Big(\frac{k_{\mathbb{43}}}{k_{\mathbb{34}}}\Big{missing})+$
$\displaystyle+$
$\displaystyle(-J^{+}_{\mathcal{C}_{3}}+-J^{-}_{\mathcal{C}_{3}}-J^{+}_{\mathcal{C}_{4}}+-J^{-}_{\mathcal{C}_{4}}-J^{+}_{\mathcal{C}_{5}}+-J^{-}_{\mathcal{C}_{5}}-J^{+}_{\mathcal{C}_{6}}+-J^{-}_{\mathcal{C}_{6}}-J^{+}_{\mathcal{C}_{7}}+-J^{-}_{\mathcal{C}_{7}}-J^{+}_{\mathcal{C}_{8}}+J^{-}_{\mathcal{C}_{8}})\ln\Big(\frac{k_{\mathbb{21}}}{k_{\mathbb{12}}}\Big{missing})+.....\Bigg{]}$
$\displaystyle=$ $\displaystyle
k_{B}\Bigg{[}J^{+}_{\mathcal{C}_{1}}\ln\Bigg(\frac{k_{\mathbb{14}}k_{\mathbb{43}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{34}}k_{\mathbb{41}}}\Bigg{missing})-J^{-}_{\mathcal{C}_{1}}\ln\Bigg(\frac{k_{\mathbb{14}}k_{\mathbb{43}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{34}}k_{\mathbb{41}}}\Bigg{missing})+.....\Bigg{]}=k_{B}\Bigg{[}J^{+}_{\mathcal{C}_{1}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}\Bigg{missing})-J^{-}_{\mathcal{C}_{1}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}\Bigg{missing})+.....\Bigg{]}$
$\displaystyle=$ $\displaystyle
k_{B}(J^{+}_{\mathcal{C}_{1}}-J^{-}_{\mathcal{C}_{1}})\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}\Bigg{missing})+.....$
$\displaystyle=$ $\displaystyle
k_{B}\sum_{\mathcal{C}}J_{\mathcal{C}}\mathcal{X}_{\mathcal{C}},\quad\text{where}\quad
J_{\mathcal{C}}=(J^{+}_{\mathcal{C}}-J^{-}_{\mathcal{C}})\quad\text{and}\quad\mathcal{X}_{\mathcal{C}}=\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}}}{\Pi^{-}_{\mathcal{C}}}\Bigg{missing}).$

The ratio of $\Pi^{\pm}_{\mathcal{C}}$ is equal to the ratio of
$J^{\pm}_{\mathcal{C}}$ for each cycle trajectory. These ratios are

$\displaystyle\frac{J^{+}_{\mathcal{C}_{1}}}{J^{-}_{\mathcal{C}_{1}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}=\frac{k_{\mathbb{43}}k_{\mathbb{31}}k_{\mathbb{14}}}{k_{\mathbb{13}}k_{\mathbb{34}}k_{\mathbb{41}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1-\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$
(121) $\displaystyle\frac{J^{+}_{\mathcal{C}_{2}}}{J^{-}_{\mathcal{C}_{2}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{2}}}{\Pi^{-}_{\mathcal{C}_{2}}}=\frac{k_{\mathbb{26}}k_{\mathbb{65}}k_{\mathbb{52}}}{k_{\mathbb{25}}k_{\mathbb{56}}k_{\mathbb{62}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1-\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$
(122) $\displaystyle\frac{J^{+}_{\mathcal{C}_{3}}}{J^{-}_{\mathcal{C}_{3}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{3}}}{\Pi^{-}_{\mathcal{C}_{3}}}=\frac{k_{\mathbb{12}}k_{\mathbb{25}}k_{\mathbb{53}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{35}}k_{\mathbb{52}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}\Big{)},$ (123)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{4}}}{J^{-}_{\mathcal{C}_{4}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{4}}}{\Pi^{-}_{\mathcal{C}_{4}}}=\frac{k_{\mathbb{12}}k_{\mathbb{26}}k_{\mathbb{65}}k_{\mathbb{53}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{35}}k_{\mathbb{56}}k_{\mathbb{62}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}-\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$ (124)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{5}}}{J^{-}_{\mathcal{C}_{5}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{5}}}{\Pi^{-}_{\mathcal{C}_{5}}}=\frac{k_{\mathbb{12}}k_{\mathbb{25}}k_{\mathbb{53}}k_{\mathbb{34}}k_{\mathbb{41}}}{k_{\mathbb{14}}k_{\mathbb{43}}k_{\mathbb{35}}k_{\mathbb{52}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}+\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$ (125)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{6}}}{J^{-}_{\mathcal{C}_{6}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{6}}}{\Pi^{-}_{\mathcal{C}_{6}}}=\frac{k_{\mathbb{12}}k_{\mathbb{26}}k_{\mathbb{64}}k_{\mathbb{41}}}{k_{\mathbb{14}}k_{\mathbb{46}}k_{\mathbb{62}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}\Big{)},$ (126)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{7}}}{J^{-}_{\mathcal{C}_{7}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{7}}}{\Pi^{-}_{\mathcal{C}_{7}}}=\frac{k_{\mathbb{12}}k_{\mathbb{25}}k_{\mathbb{56}}k_{\mathbb{64}}k_{\mathbb{41}}}{k_{\mathbb{14}}k_{\mathbb{46}}k_{\mathbb{65}}k_{\mathbb{52}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}+\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$ (127)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{8}}}{J^{-}_{\mathcal{C}_{8}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{8}}}{\Pi^{-}_{\mathcal{C}_{8}}}=\frac{k_{\mathbb{12}}k_{\mathbb{26}}k_{\mathbb{64}}k_{\mathbb{43}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{34}}k_{\mathbb{46}}k_{\mathbb{62}}k_{\mathbb{21}}}=e^{{{\rm
U}\delta T}/{k_{B}T_{0}(T_{0}+\delta T)}}e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1+\frac{{\rm U}\delta
T}{k_{B}{T_{0}}^{2}}-\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$ (128)
$\displaystyle\frac{J^{+}_{\mathcal{C}_{9}}}{J^{-}_{\mathcal{C}_{9}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{9}}}{\Pi^{-}_{\mathcal{C}_{9}}}=\frac{k_{\mathbb{34}}k_{\mathbb{46}}k_{\mathbb{65}}k_{\mathbb{53}}}{k_{\mathbb{35}}k_{\mathbb{56}}k_{\mathbb{64}}k_{\mathbb{43}}}=1,$
(129) $\displaystyle\frac{J^{+}_{\mathcal{C}_{10}}}{J^{-}_{\mathcal{C}_{10}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{10}}}{\Pi^{-}_{\mathcal{C}_{10}}}=\frac{k_{\mathbb{25}}k_{\mathbb{53}}k_{\mathbb{34}}k_{\mathbb{46}}k_{\mathbb{62}}}{k_{\mathbb{26}}k_{\mathbb{64}}k_{\mathbb{43}}k_{\mathbb{35}}k_{\mathbb{52}}}=e^{{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1+\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)},$
(130) $\displaystyle\frac{J^{+}_{\mathcal{C}_{11}}}{J^{-}_{\mathcal{C}_{11}}}$
$\displaystyle=$
$\displaystyle\frac{\Pi^{+}_{\mathcal{C}_{11}}}{\Pi^{-}_{\mathcal{C}_{11}}}=\frac{k_{\mathbb{14}}k_{\mathbb{46}}k_{\mathbb{65}}k_{\mathbb{53}}k_{\mathbb{31}}}{k_{\mathbb{13}}k_{\mathbb{35}}k_{\mathbb{56}}k_{\mathbb{64}}k_{\mathbb{41}}}=e^{-{\Delta\mu_{\rm
S}}/{k_{B}T_{0}}}\approx\Big{(}1-\frac{\Delta\mu_{\rm S}}{k_{B}T_{0}}\Big{)}.$
(131)

Using Eqs.(121)-(131) into Eq.(120), yields

$\displaystyle\dot{\sigma}$ $\displaystyle=$ $\displaystyle
J_{\mathcal{C}_{1}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{1}}}{\Pi^{-}_{\mathcal{C}_{1}}}\Bigg{missing})+J_{\mathcal{C}_{2}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{2}}}{\Pi^{-}_{\mathcal{C}_{2}}}\Bigg{missing})+J_{\mathcal{C}_{3}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{3}}}{\Pi^{-}_{\mathcal{C}_{3}}}\Bigg{missing})+J_{\mathcal{C}_{4}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{4}}}{\Pi^{-}_{\mathcal{C}_{4}}}\Bigg{missing})+J_{\mathcal{C}_{5}}\ln\Bigg(\frac{\Pi^{+}_{\mathcal{C}_{5}}}{\Pi^{-}_{\mathcal{C}_{5}}}\Bigg{missing})+.....$
(132) $\displaystyle=$ $\displaystyle
J_{\mathcal{C}_{1}}\Bigg{(}-\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}+J_{\mathcal{C}_{2}}\Bigg{(}-\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}+J_{\mathcal{C}_{3}}\Bigg{(}-\frac{{\rm U}\delta
T}{k_{B}{T_{0}}(T_{0}+\delta T)}\Bigg{)}+J_{\mathcal{C}_{4}}\Bigg{(}\frac{{\rm
U}\delta T}{k_{B}{T_{0}}(T_{0}+\delta T)}-\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}$ $\displaystyle+$ $\displaystyle
J_{\mathcal{C}_{5}}\Bigg{(}\frac{{\rm U}\delta T}{k_{B}{T_{0}}(T_{0}+\delta
T)}+\frac{\Delta\mu_{\rm
S}}{k_{B}T_{0}}\Bigg{)}+J_{\mathcal{C}_{6}}\Bigg{(}-\frac{{\rm U}\delta
T}{k_{B}{T_{0}}(T_{0}+\delta T)}\Bigg{)}+.....$ $\displaystyle=$
$\displaystyle{\rm
U}[J_{\mathcal{C}_{3}}+J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{6}}+J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}]\frac{\delta
T}{T_{0}(T_{0}+\delta
T)}+[-J_{\mathcal{C}_{1}}-J_{\mathcal{C}_{2}}-J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{7}}-J_{\mathcal{C}_{8}}+J_{\mathcal{C}_{10}}-J_{\mathcal{C}_{11}}]\frac{\Delta\mu_{\rm
S}}{T_{0}}$ $\displaystyle=$ $\displaystyle J_{\rm
E}\Big{[}\frac{1}{T_{0}}-\frac{1}{(T_{0}+\delta T)}\Big{]}+J_{\rm
S}\Big{[}\frac{\mu_{\rm L\downarrow}}{T_{0}}-\frac{\mu_{\rm
L\uparrow}}{T_{0}}\Big{]},$

where we identify the expressions of the macroscopic energy and spin currents
in terms of microscopic cycle fluxes

$\displaystyle J_{\rm E}$ $\displaystyle=$ $\displaystyle{\rm
U}[J_{\mathcal{C}_{3}}+J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{6}}+J_{\mathcal{C}_{7}}+J_{\mathcal{C}_{8}}]$
(133) $\displaystyle J_{\rm S}$ $\displaystyle=$
$\displaystyle[-J_{\mathcal{C}_{1}}-J_{\mathcal{C}_{2}}-J_{\mathcal{C}_{4}}+J_{\mathcal{C}_{5}}+J_{\mathcal{C}_{7}}-J_{\mathcal{C}_{8}}+J_{\mathcal{C}_{10}}-J_{\mathcal{C}_{11}}],$
(134)

which are equivalent to Eqs. (62) and (60) of the main text.