††thanks<EMAIL_ADDRESS>
# Thermodynamic theory of inverse current in coupled quantum transport

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

###### Abstract

The inverse current in coupled (ICC) quantum transport, where one induced
current opposes all thermodynamic forces of a system, is a highly counter-
intuitive transport phenomenon. Using an exactly solvable model of strongly-
coupled quantum dots, we present thermodynamic description of ICC in energy
and spin-induced particle currents, with potential applications towards
unconventional and autonomous nanoscale thermoelectric generators. Our
analysis reveals the connection between microscopic and macroscopic
formulations of entropy production rates, elucidating the often-overlooked
role of proper thermodynamic forces and conjugate fluxes in characterizing
genuine ICC. In our model, the seemingly paradoxical results of ICC in the
energy current arise from chemical work done by current-carrying quantum
particles, while in spin-induced particle current, it stems from the relative
competition between electron reservoirs controlling one particular transition.

##  I Introduction

Force and flux are two central thermodynamic quantities interlinked through a
cause-and-effect relationship. Thermodynamic force emerges due to a gradient
in the system and engenders a response in the form of flux within the system,
causing a departure from equilibrium. Consequently, the resultant flux aligns
in the direction of the applied force to establish a new equilibrium. One of
the most intriguing transport phenomena is absolute negative mobility (ANM),
wherein the system’s response i.e., the generated current, operates against
the driving force Cleuren and den Broeck (2001); Cleuren and Van den Broeck
(2002); Eichhorn  _et al._ (2002). The absence of ANM, at near-equilibrium for
a single flux-force pair, can be understood from the sign of the entropy
production rate ($\dot{\Sigma}=J\mathcal{F}$) to remain consistent with the
second law, where $J$ represents the thermodynamic flux conjugate to the
thermodynamic force $\mathcal{F}$. Initially, ANM was considered as a
consequence of the quantum effect Höpfel  _et al._ (1986); Keay  _et al._
(1995). However, theoretical works revealing the lack of any fundamental laws
prohibiting ANM in systems far from equilibrium created significant interests
Cleuren and den Broeck (2001); Cleuren and Van den Broeck (2002); Eichhorn
_et al._ (2002); Cleuren and Van den Broeck (2003); Machura  _et al._ (2007);
Dandogbessi and Kenfack (2015); Cividini  _et al._ (2018); Słapik  _et al._
(2019), which have been verified in experiments Keay  _et al._ (1995); Ros
_et al._ (2005); Nagel  _et al._ (2008).

The above situation changes dramatically in the presence of coupled transports
with multiple force-flux pairs, where the entropy production rate near
equilibrium takes the general form
$\dot{\Sigma}=\sum_{i}J_{i}\mathcal{F}_{i}$. This could in principle lead to
the concept of Inverse Current in Coupled transport (ICC), where a current
flows against all forces, including its own conjugate force. It is crucial to
distinguish ICC from a cross-effect, where a current is produced in response
to a non-conjugate force, despite the fact that its own conjugate force is set
to zero. While a cross-effect (Seebeck and Peltier) is a regular outcome in
coupled transport, ICC involves a current opposing to all forces present in
the system, including its own conjugate force, appearing highly
counterintuitive, yet not forbidden as long as the overall entropy production
rate remains in the positive domain. To examine ICC, one must therefore
consider at least two nonzero forces and their corresponding fluxes. For
example, the thermal (energy) force $\mathcal{F}_{\rm{E}}$, resulting from a
temperature gradient ($\Delta T$) and the particle force
$\mathcal{F}_{\rm{N}}$, arising from chemical potential gradient
($\Delta\mu$), could drive both particle $(J_{\rm{N}})$ and energy
$(J_{\rm{E}})$ currents. As a result, entropy production rate is expressed as
$\dot{\Sigma}=J_{\rm{E}}\mathcal{F}_{\rm{E}}+J_{\rm{N}}\mathcal{F}_{\rm{N}}$.
When both currents may counter-operate their non-conjugate forces, for ICC to
occur, either $J_{\rm{E}}$ or $J_{\rm{N}}$ must exhibit a sign opposite to
both forces. The immediate consequence is that simultaneous ICC in both fluxes
is impossible, as it would violate the second law of thermodynamics, resulting
in a negative entropy production rate.

Recently, based on a classical Lieb-Liniger model Lieb and Liniger (1963),
characterized by 1D-interacting Hamiltonian of a diatomic gas, a potential
platform for ICC, separately for energy and particle currents, has been
demonstrated by Wang et al. Wang  _et al._ (2020). The seemingly paradoxical
result in the above classical model originates from the occurrence of self-
organization within the system in response to applied forces. Since, the
general criteria for ICC are yet unknown either in classical/quantum settings,
showcasing the existence of ICC in quantum systems, still remains a difficult
task. Only a few attempts are made Zhang and Xie (2021); Zhang  _et al._
(2023) which however fell short of establishing ICC in terms of proper
thermodynamic forces. Though numerical results managed to produce the flow of
one current against both gradients ($\Delta T$ and $\Delta\mu$), these
observations do not qualify as ICC, as thermodynamic forces and gradients are
altogether distinct quantities. Thus, realizing genuine ICC in quantum systems
with appropriate thermodynamic analysis still remains a highly challenging
endeavor.

In this paper, we provide a resolution to this problem by constructing a
thermodynamic theory of inverse current using an exactly solvable model of a
three-terminal Coulomb-coupled quantum dots (QDs), a simple variant of the
extensively studied Sánchez-Büttiker model Sánchez and Büttiker (2011),
explored in number of contexts, such as quantum transport Timm (2007);
Walschaers  _et al._ (2013); Kulvelis  _et al._ (2015); Thingna  _et al._
(2016); Maier  _et al._ (2019); Thingna  _et al._ (2020); Ehrhardt  _et al._
(2021), quantum information Strasberg  _et al._ (2013); Kutvonen  _et al._
(2016); Ptaszyński and Esposito (2019) and thermoelectricity Esposito  _et
al._ (2009); Nakpathomkun  _et al._ (2010); Sothmann and Büttiker (2012);
Donsa  _et al._ (2014); Zhang  _et al._ (2015); Thierschmann  _et al._ (2015);
Sierra  _et al._ (2016); Whitney  _et al._ (2016); Erdman  _et al._ (2017,
2018); Whitney  _et al._ (2018); Mazal  _et al._ (2019); Manzano  _et al._
(2020); Wang  _et al._ (2022), thermal rectification Ruokola and Ojanen
(2011); Werlang  _et al._ (2014); Hartmann  _et al._ (2015); Zhang  _et al._
(2017); Aligia  _et al._ (2020); Tesser  _et al._ (2022); Ghosh  _et al._
(2022), and many others Esposito  _et al._ (2012); Jiang  _et al._ (2015);
Yang  _et al._ (2019). Our findings unravel the crucial role played by the QD
interaction in establishing ICC in energy and particle currents. The key
element of our analysis lies in identifying macroscopic and microscopic
formulas for thermodynamic forces and fluxes. Most importantly, we have
established conditions for genuine ICC in both energy and (spin-polarized)
particle currents, without violating the second law of thermodynamics.

The structure of the paper is as follows: Section II outlines our model, which
is based on Coulomb-coupled QDs. We explore the system’s dynamics in Sec. III,
analyze steady-state currents in Sec. IV, and examine entropy production,
including macroscopic and microscopic forces and fluxes, in Sec. V. Section VI
provides a summary of our results. Finally, we conclude in Sec. VII.

##  II Theoretical Model

Our model consists of two strongly and capacitively coupled quantum dots
(${\rm QDs}$), labeled by left (${\rm QD_{L}}$) and right (${\rm QD_{R}}$),
respectively [FIG. 1]. These dots solely interact via a long-range Coulomb
force restricting particle exchange due to Coulomb blockade Ruokola and Ojanen
(2011); Zhang  _et al._ (2017); Ghosh  _et al._ (2022) while facilitating
energy exchange within the dots through Coulomb interaction
($\kappa_{\rm{c}}$). The ${\rm QD_{L}}$ is simultaneously tunnel-coupled with
two fermionic reservoirs Sánchez and Büttiker (2011); Ghosh  _et al._ (2012);
Zhang  _et al._ (2018); Ghosh  _et al._ (2022), labeled above (a) and below
(b), allowing particle flow between the two reservoirs. The ${\rm QD_{R}}$ is
tunnel-coupled with only one reservoir (r), permitting both particle and
energy exchange between the dot and the coupled lead. While for spinless
electrons, $\kappa_{\rm{c}}$ is always positive, w.l.o.g, we consider spin-
polarized Agrawal and Solomon (2004); Vandaele  _et al._ (2017); Gupt  _et
al._ (2021) electrons to account for a more comprehensive study. To facilitate
electron exchange between the leads a and b through ${\rm QD_{L}}$, we assume
that reservoirs a and b consist of only spin-down ($\downarrow$) electrons,
while reservoir r is filled with spin-up ($\uparrow$) electrons. The setup
could lead to both attractive (negative) and repulsive (positive) interaction
between the ${\rm QDs}$ depending on the relative strength of the $\kappa_{\rm
c}$ and the spin-spin interaction ($\kappa_{\rm s}$). The Hamiltonian
governing the coupled ${\rm QD}$ system is given by Werlang  _et al._ (2014);
Ghosh  _et al._ (2022); Gupt  _et al._ (2023, 2024),

$\begin{split}H_{\rm{s}}=\varepsilon_{\rm{L}}\mathcal{N}_{\rm{L}\downarrow}+\varepsilon_{\rm{R}}\mathcal{N}_{\rm{R}\uparrow}+\kappa_{c}\mathcal{N}_{\rm{L}\downarrow}\mathcal{N}_{\rm{R}\uparrow}+\kappa_{s}\sigma^{z}_{\rm{L}\downarrow}\sigma^{z}_{\rm{R}\uparrow}.\end{split}$
(1)

In the above Hamiltonian, $\varepsilon_{\alpha}$ $(\alpha=\rm{L},\rm{R})$,
represents the single-particle energy level associated with the $\alpha$’th
${\rm QD}$. Due to the Coulomb blockade, the electron density in the dots is
low, restricting occupancy to either zero or one. Under this condition, the
eigenstates of ${\rm QD_{\alpha}}$ are either $|0\rangle$ or
$\ket{\downarrow}(\ket{\uparrow})$, with energy eigenvalues $0$ and
$\varepsilon_{\rm{L}}(\varepsilon_{\rm{R}})$, respectively. The corresponding
number operators for ${\rm QD_{\alpha}}$ are
$\mathcal{N}_{\rm{L}\downarrow}=d^{\dagger}_{\rm{L}\downarrow}d_{\rm{L}\downarrow}$
and
$\mathcal{N}_{\rm{R}\uparrow}=d^{\dagger}_{\rm{R}\uparrow}d_{\rm{R}\uparrow}$,
where $d^{\dagger}_{\rm{L}\downarrow}$ ($d_{\rm{L}\downarrow}$) and
$d^{\dagger}_{{\rm{R}}\uparrow}$ ($d_{{\rm{R}}\uparrow}$) denote the electron
creation (annihilation) operators for the respective ${\rm QD_{\alpha}}$,
obeying to the anti-commutation relations
$\\{d_{\rm{L}\downarrow},d^{\dagger}_{\rm{L}\downarrow}\\}=1=\\{d_{{\rm{R}}\uparrow},d^{\dagger}_{{\rm{R}}\uparrow}\\}$.
The last term in Eq. (1) represents the spin-spin interaction energy where
$\sigma^{z}_{\rm{L\downarrow(R\uparrow)}}=1\mp
2\mathcal{N}_{\rm{L\downarrow(R\uparrow)}}$, satisfying standard operator
algebra $\sigma^{z}_{\rm{L}\downarrow}|\downarrow\rangle=-1|\downarrow\rangle$
and $\sigma^{z}_{\rm{R}\uparrow}|\downarrow\rangle=+1|\uparrow\rangle$.

As the two QDs are strongly and capacitively coupled, the overall system
Hamiltonian is diagonal in the eigenbasis of the individual ${\rm QD}$. This
can be represented by the tensor product of the number operator’s eigenbasis
of the coupled ${\rm QDs}$. For convenience, the four eigenstates
$\\{|0\rangle,\ket{\downarrow}\\}\otimes\\{|0\rangle,\ket{\uparrow}\\}$, are
labeled by $|\mathbb{A}\rangle=\ket{00}$, $|\mathbb{B}\rangle=\ket{\downarrow
0}$, $|\mathbb{C}\rangle=\ket{0\uparrow}$,
$|\mathbb{D}\rangle=\ket{\downarrow\uparrow}$ and their corresponding
eigenenergies ($\varepsilon_{\mathbb{i}}$, $\mathbb{i=A,B,C,D}$) are
$\varepsilon_{\mathbb{A}}=0$, $\varepsilon_{\mathbb{B}}=\varepsilon_{\rm{L}}$,
$\varepsilon_{\mathbb{C}}=\varepsilon_{\rm{R}}$ and
$\varepsilon_{\mathbb{D}}=\varepsilon_{\rm{L}}+\varepsilon_{\rm{R}}+\kappa$
respectively [FIG. 1: Inset]. Consequently, the energy of the most excited
state of the composite system is
$\varepsilon_{\rm{L}}+\varepsilon_{\rm{R}}+\kappa$, where, we assume w.l.o.g,
$\varepsilon_{\rm{L}}<\varepsilon_{\rm{R}}$. Although, $\kappa_{\rm{c}}$ and
$\kappa_{\rm{s}}$ are non-negative quantities, the overall interaction
$\kappa=\kappa_{\rm{c}}-\kappa_{\rm{s}}$, can take any real values involving
positive and negative domains of interaction. It is important to note that
when $\kappa$ is negative and $|\kappa|>\varepsilon_{\rm L}$, the energy
states $|\mathbb{C}\rangle$ and $|\mathbb{D}\rangle$ swap their positions
[FIG. 1: Inset]. This will play a crucial role in exhibiting ICC behavior, as
explored in Sec. VI.

Figure 1: Main: (a) Schematic diagram of the three-terminal Coulomb-coupled
QDs. Inset: The energy level diagram of the eigenstates of coupled the QDs for
$\kappa>0$ (left) and $\kappa<0;|\kappa|>\varepsilon_{\rm{L}}$ (right). From
the general model (a), we can construct two reduced models:- (b) Reduced
model-I, by considering the temperatures of the reservoirs $(a)$ and $(b)$
equal. (c) Reduced model-II, by setting up the temperatures of the reservoirs
$(a)$ and $(r)$ the same.

Finally, the reservoirs are dense with fermionic particles, i.e., electrons,
characterized by temperature and chemical potential. The Hamiltonian of the
$\lambda$’th spin-polarized reservoir with electron spin
$\sigma=\\{\uparrow,\downarrow\\}$, is defined as
$H_{\rm{B}}^{\lambda\sigma}=\sum_{k}(\epsilon^{\lambda\sigma}_{k}-\mu_{\lambda\sigma})c_{{\lambda\sigma}k}^{\dagger}c_{{\lambda\sigma}k}$,
where $\epsilon^{\lambda}_{k}$ is the energy of the non-interacting electrons
for the reservoir $\lambda$, $k$ is the continuous wavenumber,
$\mu_{\lambda\sigma}$ is the chemical potential, and
$c^{\dagger}_{\lambda\sigma}(c_{\lambda\sigma})$ represents the creation
(annihilation) operator. The total Hamiltonian of all three reservoirs is
given by
$H_{\rm{B}}=H_{\rm{B}}^{a\downarrow}+H_{\rm{B}}^{b\downarrow}+H_{\rm{B}}^{r\uparrow}$.
The ${\rm QD_{R}}$ (${\rm QD_{L}}$) is weakly coupled to the reservoir(s) $r$
($a$ and $b$) to ensure sequential tunneling Zhang  _et al._ (2017, 2018);
Ghosh  _et al._ (2022), such that only one ${\rm QD}$ at a time is involved
with particle tunneling with the coupled lead. The tunnel-coupled Hamiltonians
are characterized by the coupling constant $t_{k}^{\alpha\sigma\lambda}$ and
are given by

$\displaystyle H_{\rm{T}}^{{\rm{L}\downarrow}{a(b)}}$ $\displaystyle=$
$\displaystyle\hbar\sum_{k}[t^{{\rm{L}\downarrow}{a(b)}}_{k}c_{a(b)\downarrow
k}^{\dagger}d_{{\rm{L}}\downarrow}+t^{{\rm{L}\downarrow}{a(b)*}}_{k}d_{{\rm{L}}\downarrow}^{\dagger}c_{a(b)\downarrow
k}],$ $\displaystyle H_{{\rm{T}}}^{{\rm{R}\uparrow}{r}}$ $\displaystyle=$
$\displaystyle\hbar\sum_{k}[t^{{\rm{R}\uparrow}{r}}_{k}c_{r\uparrow
k}^{\dagger}d_{{\rm{R}}\uparrow}+t^{{\rm{R}\uparrow}{r*}}_{k}d_{{\rm{R}}\uparrow}^{\dagger}c_{r\uparrow
k}].$ (2)

where, $H_{\rm{T}}^{{\rm{L}\downarrow}{a(b)}}$ signifies the interaction
between $\rm{QD_{L}}$ and reservoir $a(b)$ with the exchange of $\downarrow$
electrons and $H_{\rm{T}}^{{\rm{R}\uparrow}{r}}$ represents the interaction
between $\rm{QD_{R}}$ and reservoir $r$ through the exchange of $\uparrow$
electrons. Further, the transitions between
$|\mathbb{B}\rangle\leftrightarrow|\mathbb{C}\rangle$ and
$|\mathbb{A}\rangle\leftrightarrow|\mathbb{D}\rangle$ are prohibited due to
Coulomb blockade and the sequential tunneling approximation Zhang  _et al._
(2017, 2018); Ghosh  _et al._ (2022), respectively. So, there are in total
four authorized transitions between the eigenstates of the composite system,
and each of the three reservoirs $a$, $b$, and $r$ controls two transitions,
i.e., $|\mathbb{A}\rangle\leftrightarrow|\mathbb{B}\rangle$ and
$|\mathbb{C}\rangle\leftrightarrow|\mathbb{D}\rangle$ are triggered by leads
$a$ and/or $b$, while reservoir $r$ governs the transitions between
$|\mathbb{A}\rangle\leftrightarrow|\mathbb{C}\rangle$ and
$|\mathbb{B}\rangle\leftrightarrow|\mathbb{D}\rangle$. The transition energy
between the eigenstates $|\mathbb{i}\rangle\rightarrow|\mathbb{j}\rangle$ is
defined as
$\omega_{\mathbb{ij}}=\varepsilon_{\mathbb{j}}-\varepsilon_{\mathbb{i}}$,
where $\varepsilon_{\mathbb{i}}$ is the eigenenergy of $H_{\rm{s}}$ with the
eigenstate $\\{|\mathbb{i}\rangle\\}$ and the four transition energies are
respectively $\omega_{\mathbb{AB}}=\varepsilon_{\rm{L}}$,
$\omega_{\mathbb{AC}}=\varepsilon_{\rm{R}}$,
$\omega_{\mathbb{CD}}=\varepsilon_{\rm{L}}+\kappa$, and
$\omega_{\mathbb{BD}}=\varepsilon_{\rm{R}}+\kappa$. To investigate the
system’s dynamics, we will next evaluate the transition rates using the
Lindblad Master Equation (LME) under the Born-Markov approximation.

##  III Evaluation of rate equations using LME

The state of the composite system is described by the reduced density matrix
$\rho_{\rm{s}}(t)=\Tr_{\rm{B}}\\{\rho_{\rm{tot}}(t)\\}$, where
$\rho_{\rm{tot}}(t)$ is the total density matrix of the system and bath
combined. We utilize strong-coupling formalism to derive the LME governing the
time evolution of the reduced density matrix $\rho_{\rm{s}}(t)$ of the system,
under the Born, Markov, and Secular (BMS) approximation Breuer and Petruccione
(2007); Ghosh  _et al._ (2022); Strasberg (2022); Gupt  _et al._ (2023)
(Appendix A):

$\frac{d}{dt}\rho_{\rm{s}}(t)=\sum_{\lambda}\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)];\quad\lambda=a,b,r.$
(3)

It is worth noting that the strong coupling formalism pertains to the
interaction between the QDs, while we maintain the assumption of weak coupling
between the system and its surrounding environment. This allows us to safely
implement the BMS approximation and derive the above LME [Eq. (3)] purely
based on the eigenstates of the full system Hamiltonian $H_{\rm{s}}$. As a
result, the dissipation mechanism of each QD is influenced not only by its
coupling to its individual bath but also by the interactions between QDs
themselves. This consideration appears crucial for accurately characterizing
heat flow and obtaining outcomes across a wide range of system parameters, as
examined below. $\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]$ in the above
equation is the Lindblad super-operator, defined as

$\begin{split}\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]=&\sum_{\\{\omega_{\alpha}\\}>0}\left\\{\gamma_{\lambda}(\omega_{\alpha})f^{+}_{\lambda}(\omega_{\alpha})\left[d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})\rho_{\rm{s}}d_{\alpha\sigma}(\omega_{\alpha})-\frac{1}{2}{\\{\rho_{\rm{s}}},d_{\alpha\sigma}(\omega_{\alpha})d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})\\}\right]\right.\\\
&\left.+\gamma_{\lambda}(\omega_{\alpha}){f}^{-}_{\lambda}(\omega_{\alpha})\left[d_{\alpha\sigma}(\omega_{\alpha})\rho_{\rm{s}}d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})-\frac{1}{2}{\\{\rho_{\rm{s}}},d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})d_{\alpha\sigma}(\omega_{\alpha})\\}\right]\right\\},\end{split}$
(4)

where $\gamma_{\lambda}(\omega_{\alpha})$ is the bare electron transfer rate
between the reservoir $\lambda$ and coupled ${\rm QD_{\alpha}}$. The explicit
forms in terms of the system-reservoir coupling constants can be calculated
using Fermi’s golden rule, as
$\gamma_{\lambda}(\omega_{\alpha})\equiv\gamma_{\lambda\sigma}(\omega_{\alpha})=2\pi\sum_{k}|t_{k}^{\alpha\sigma\lambda}|^{2}\delta\big{(}\omega_{\alpha}-\epsilon_{k}^{\lambda\sigma}\big{)}$,
where $\omega_{\alpha}$ represents the required amount of energy associated
with the transition between ${\rm QD_{\alpha}}$ and its coupled lead. For
simplicity of notations, we use
$\gamma_{a|b}(\omega_{\rm{L}})\equiv\gamma_{a|b}$ and
$\gamma_{r}(\omega_{\rm{R}})\equiv\gamma_{r}$, where
$\omega_{\rm{L}}=\\{\omega_{\mathbb{AB}},\omega_{\mathbb{CD}}\\}$ and
$\omega_{\rm{R}}=\\{\omega_{\mathbb{AC}},\omega_{\mathbb{BD}}\\}$, if the
particle is entering (excitation) into the system, and similarly, for particle
de-excitation (leaving) w.r.t. system,
$\omega_{\rm{L}}=\\{\omega_{\mathbb{BA}},\omega_{\mathbb{DC}}\\}$ and
$\omega_{\rm{R}}=\\{\omega_{\mathbb{CA}},\omega_{\mathbb{DB}}\\}$. Here, the
symbol $a|b$ denotes the transition controlled by either $a$ or $b$ bath.
Finally, $f^{\pm}_{\lambda}(\omega_{\mathbb{ij}})$ represents the Fermi
distribution function (FDF) related to the transition
$|\mathbb{i}\rangle\rightarrow|\mathbb{j}\rangle$ driven by the $\lambda$’th
reservoir which cost $\omega_{\mathbb{ij}}$ amount of energy. The notation
‘$\pm$’ which is used throughout the text Sánchez and Büttiker (2011); Zhang
and Xie (2021), where the ‘$\pm$’ sign refers to particle excitation (entering
into the system) and the ‘$-$’ sign refers to particle de-excitation (leaving
from the system) w.r.t the system. Hence, if the transition from
$|\mathbb{i}\rangle$ to $|\mathbb{j}\rangle$ is governed by the particle
excitation w.r.t the system then the reverse process would signify the
tunneling of the particle from the system to the environment. Thus, we
establish the following relation between the FDFs associated with the above
transitions Ghosh  _et al._ (2022); Tesser  _et al._ (2022),

$f^{+}_{\lambda}(\omega_{\mathbb{ij}})=1-f^{-}_{\lambda}(\omega_{\mathbb{ji}})=\bigg{[}{1+\exp\bigg({\frac{\omega_{\mathbb{ij}}-\mu_{\lambda}}{k_{B}T_{\lambda}}})}\bigg{]}^{-1},$
(5)

where $T_{\lambda}$ and $\mu_{\lambda}$ are respectively the temperature and
chemical potential of the fermionic reservoir $\lambda$. For the sake of
convenience, in Eq. (5), we have dropped the $\sigma$ from the temperature and
chemical potential. Anticipating what will become clear as we proceed, the
above notation has the advantage that it does not rely on the specific energy
level ordering of the system eigenstates. Therefore, all subsequent
calculations can be done in a compact form without referring to any particular
arrangement of the system eigen-energies.

So, the time evolution of the occupation probabilities which are the diagonal
elements of the reduced density matrix, i.e.,
$\rho_{\mathbb{ii}}=\langle\mathbb{i}|\rho_{\rm{s}}(t)|\mathbb{i}\rangle$, can
be evaluated from the LME [Eq. (3)] as

$\begin{split}\dot{\rho}_{\mathbb{i}}\equiv\frac{d\rho_{\mathbb{ii}}}{dt}=\langle\mathbb{i}|\dot{\rho_{\rm{s}}}(t)|\mathbb{i}\rangle=\sum_{\lambda=a,b,r}\langle\mathbb{i}|\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]|\mathbb{i}\rangle.\end{split}$
(6)

Using the Lindbladians defined in Eq. (4), we can determine the rate equations
as follows

$\begin{split}\dot{\rho}_{\mathbb{A}}=&-\Gamma^{a+}_{\mathbb{AB}}-\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{r+}_{\mathbb{AC}}\equiv-\Gamma^{ab+}_{\mathbb{AB}}-\Gamma^{r+}_{\mathbb{AC}};\\\
\dot{\rho}_{\mathbb{B}}=&\Gamma^{a+}_{\mathbb{AB}}+\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{r+}_{\mathbb{BD}}\equiv\Gamma^{ab+}_{\mathbb{AB}}-\Gamma^{r+}_{\mathbb{BD}};\\\
\dot{\rho}_{\mathbb{C}}=&-\Gamma^{a+}_{\mathbb{CD}}-\Gamma^{b+}_{\mathbb{CD}}+\Gamma^{r+}_{\mathbb{AC}}\equiv-\Gamma^{ab+}_{\mathbb{CD}}+\Gamma^{r+}_{\mathbb{AC}};\\\
\dot{\rho}_{\mathbb{D}}=&\Gamma^{a+}_{\mathbb{CD}}+\Gamma^{b+}_{\mathbb{CD}}+\Gamma^{r+}_{\mathbb{BD}}\equiv\Gamma^{ab+}_{\mathbb{CD}}+\Gamma^{r+}_{\mathbb{BD}},\end{split}$
(7)

where,
$\Gamma^{ab\pm}_{\mathbb{ij}}\equiv\Gamma^{a\pm}_{\mathbb{ij}}+\Gamma^{b\pm}_{\mathbb{ij}}$.
The net transition rate $\Gamma^{\lambda\pm}_{\mathbb{ij}}$ from level
$|\mathbb{i}\rangle$ to $|\mathbb{j}\rangle$ mediated by the reservoir
$\lambda$ ($\lambda=a,b,r$) for both the ‘$\pm$’ processes can be expressed in
a compact form as

$\begin{split}\Gamma_{\mathbb{ij}}^{\lambda\pm}\equiv\Gamma_{\mathbb{i\mapsto
j}}^{\lambda\pm}=&\gamma_{\lambda}f^{\pm}_{\lambda}(\omega_{\mathbb{ij}})\rho_{\mathbb{i}}-\gamma_{\lambda}f^{\mp}_{\lambda}(\omega_{\mathbb{ji}})\rho_{\mathbb{j}},\\\
=&{\rm{k}}_{\mathbb{ij}}^{\lambda\pm}\rho_{\mathbb{i}}-{\rm{k}}_{\mathbb{ji}}^{\lambda\mp}\rho_{\mathbb{j}}\end{split}$
(8)

where $\Gamma_{\mathbb{ij}}^{\lambda+}=-\Gamma_{\mathbb{ji}}^{\lambda-}$ and

${{\rm{k}}}_{\mathbb{ij}}^{\lambda\pm}=\gamma_{\lambda}f^{\pm}_{\lambda}(\omega_{\mathbb{ij}})\quad;\quad{{\rm{k}}}_{\mathbb{ji}}^{\lambda\mp}=\gamma_{\lambda}f^{\mp}_{\lambda}(\omega_{\mathbb{ji}}).$
(9)

Using Eq. (5), one can verify that the above rate coefficients for fermionic
reservoirs satisfy an interesting relation

${{\rm{k}}}_{\mathbb{ij}}^{\lambda\pm}+{{\rm{k}}}_{\mathbb{ji}}^{\lambda\mp}=\gamma_{\lambda}.$
(10)

For simplicity, we further assume all $\gamma_{\lambda}$ are equal when
calculating the closed-form analytical expression for the steady-state
transition rates in Sec. IV.

##  IV Steady State Currents and Dynamics

In this section, we derive the expressions for steady-state currents under the
grand canonical formalism, since the composite system weakly interacts with
multiple reservoirs and exchanges both the energy and particles. The initial
density operator of the system $\rho_{\rm{s}}(0)$ at equilibrium can be
expressed as Strasberg (2022)

$\rho_{\rm{s}}^{\rm{eq}}=\frac{e^{-\beta(H_{\rm{s}}-\mu\mathcal{N})}}{\mathcal{Z(\beta,\mu)}},$
(11)

where $\mathcal{Z(\beta,\mu)}=\Tr[e^{-\beta(H_{\rm{s}}-\mu\mathcal{N})}]$ is
the grand canonical partition function and
$\mathcal{N}=\mathcal{N}_{\rm{L}\downarrow}+\mathcal{N}_{\rm{R}\uparrow}$
being the total particle number operator of the two QDs. The $\beta$ and $\mu$
are the effective inverse temperature and chemical potential,
$\beta=\sum_{\lambda}\beta_{\lambda}$ and $\mu=\sum_{\lambda}\mu_{\lambda}$
respectively. We assume that environmental interaction slightly perturbs the
system from its initial equilibrium state
$\rho_{\rm{s}}(0)=\rho_{\rm{s}}^{\rm{eq}}$ to $\rho_{\rm{s}}(t)$, such that
$\delta\rho_{\rm{s}}(t)=\rho_{\rm{s}}(t)-\rho_{\rm{s}}(0)\equiv\mathcal{O}(\xi)$,
where $\xi$ is a small expansion parameter. Equating system von-Neumann
entropy (times $k_{B}$) with thermodynamic entropy near close to equilibrium
Schnakenberg (1976); Esposito  _et al._ (2010); Landi and Paternostro (2021)

$\begin{split}\mathcal{S}_{\rm{s}}(t)=-k_{B}\Tr_{\rm{s}}[\rho_{\rm{s}}(t)\ln{\rho_{\rm{s}}(t)}],\\\
\end{split}$ (12)

one can derive the following expression for
$\varDelta\mathcal{S}_{\rm{s}}(t)=\mathcal{S}_{\rm{s}}(t)-\mathcal{S}_{\rm{s}}(0)$,
considering up to first order change in $\xi$,

$\begin{split}\varDelta\mathcal{S}_{\rm{s}}(t)=&-k_{B}\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)\ln{\rho_{\rm{s}}(0)}]-k_{B}\Tr_{\rm{s}}\bigg{[}\delta\rho_{\rm{s}}(t)\ln\left\\{1+\frac{\delta\rho_{\rm{s}}(t)}{\rho_{\rm{s}}(0)}\right\\}\bigg{]}-k_{B}\Tr_{\rm{s}}\bigg{[}\rho_{\rm{s}}(0)\ln\left\\{1+\frac{\delta\rho_{\rm{s}}(t)}{\rho_{\rm{s}}(0)}\right\\}\bigg{]}\\\
=&-k_{B}\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)\ln{\rho_{\rm{s}}^{\rm
eq}}]+\mathcal{O}(\xi^{2})\\\
\equiv&\;k_{B}\beta\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)H_{\rm{s}}]-k_{B}\beta\mu\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)\mathcal{N}].\end{split}$
(13)

Figure 2: Schematic representation of the clockwise (Left) and anti-clockwise
(Right) transition cycles, induced by energy flow from one bath to another,
associated with the particle exchange between QDs and the coupled reservoirs.
At steady state, all transition rates of the clockwise (anti-clockwise) cycles
are equal to each other and represented by
$\Gamma_{\circlearrowright}(\Gamma_{\circlearrowleft})$ [Cf. Eq. (21)].

Comparing the above equation and the first law of thermodynamics in the
presence of chemical work done, we identify the change in energy and particle
number as follows

$\varDelta
E=\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)H_{\rm{s}}]\quad;\quad\varDelta
N=\Tr_{\rm{s}}[\delta\rho_{\rm{s}}(t)\mathcal{N}].$ (14)

Noting $\delta\rho_{\rm{s}}(t)=\rho_{\rm{s}}(t)-\rho_{\rm{s}}(0)$, and
$H_{\rm{s}}$, $\mathcal{N}$ are time-independent, we immediately recognize the
net energy flux ($J_{\rm{E}}$) and particle flux ($J_{\rm{N}}$) as

$\begin{split}J_{\rm{E}}(t)=&\Tr_{\rm{s}}[\dot{\rho}_{\rm{s}}(t)H_{\rm{s}}]=\sum_{\lambda=a,b,r}\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]H_{\rm{s}}],\\\
J_{\rm{N}}(t)=&\Tr_{\rm{s}}[\dot{\rho}_{\rm{s}}(t)\mathcal{N}]=\sum_{\lambda=a,b,r}\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]\mathcal{N}],\end{split}$
(15)

where we use Eq. (3) for $\dot{\rho}_{\rm{s}}(t)$. Identifying $J_{\rm{E(N)}}$
as the sum of the contributions of energy (particle) flux associated with all
three reservoirs, we can write

$\displaystyle
J_{\rm{E}}(t)=\sum_{\lambda=a,b,r}J^{\lambda}_{\rm{E}}(t)\quad;\quad
J_{\rm{N}}(t)=\sum_{\lambda=a,b,r}J^{\lambda}_{\rm{N}}(t),$ (16)

where $J^{\lambda}_{\rm{E(N)}}$ is positive if energy (particle) current flows
from the reservoir $\lambda$ to the system. By comparing Eq.(15) and Eq.(16),
we thus find

$J^{\lambda}_{\rm{E}}(t)=\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho(t)]H_{\rm{s}}]\;;\;J^{\lambda}_{\rm{N}}(t)=\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho(t)]\mathcal{N}].$
(17)

Similarly, we can evaluate the expression of heat current due to the
$\lambda$’th reservoir as

$J^{\lambda}_{\rm{Q}}(t)=J^{\lambda}_{\rm{E}}(t)-\mu_{\lambda}J^{\lambda}_{\rm{N}}(t).$
(18)

Now, at the steady state, both energy and particle currents would be
independent of the time, i.e., Eq. (18) reduces to

$J^{\lambda}_{\rm{Q}}=J^{\lambda}_{\rm{E}}-\mu_{\lambda}J^{\lambda}_{\rm{N}}=\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho_{ss}]H_{\rm{s}}]-\mu_{\lambda}\Tr_{\rm{s}}[\mathcal{L}_{\lambda}[\rho_{ss}]\mathcal{N}].$
(19)

Using Eq. (4) we obtain the form of $\mathcal{L}_{\lambda}[\rho_{ss}]$, which
leads to the formal expressions for the energy, particle, and heat currents as
follows Kosloff (2013)

$\begin{split}J^{\lambda}_{\rm{E}}=&\sum_{\\{\omega_{\mathbb{ij}}\\}}{\omega_{\mathbb{ij}}}\Gamma_{\mathbb{ij}}^{\lambda+}=\sum_{\\{\omega_{\mathbb{ji}}\\}}{\omega_{\mathbb{ji}}}\Gamma_{\mathbb{ji}}^{\lambda-};\\\
J^{\lambda}_{\rm{N}}=&\sum_{\\{\omega_{\mathbb{ij}}\\}}\Gamma_{\mathbb{ij}}^{\lambda+}=\sum_{\\{\omega_{\mathbb{ji}}\\}}\Gamma_{\mathbb{ji}}^{\lambda-};\\\
J^{\lambda}_{\rm{Q}}=&\sum_{\\{\omega_{\mathbb{ij}}\\}}({\omega_{\mathbb{ij}}}-\mu_{\lambda})\Gamma_{\mathbb{ij}}^{\lambda+}=\sum_{\\{\omega_{\mathbb{ji}}\\}}({\omega_{\mathbb{ji}}}-\mu_{\lambda})\Gamma_{\mathbb{ji}}^{\lambda-}.\end{split}$
(20)

The explicit expressions for all three currents are given in Appendix B. At
the steady state, the populations of the various eigenstates become time-
invariant, i.e.,
$\dot{\rho}_{\mathbb{A}}=\dot{\rho}_{\mathbb{B}}=\dot{\rho}_{\mathbb{C}}=\dot{\rho}_{\mathbb{D}}=0$;
therefore, we obtain

$\begin{split}\Gamma^{ab+}_{\mathbb{AB}}=\Gamma^{r+}_{\mathbb{BD}}=\Gamma^{ab-}_{\mathbb{DC}}=\Gamma^{r-}_{\mathbb{CA}}=\Gamma_{\mathbb{ABCDA}}\equiv\Gamma_{\circlearrowright};\\\
\Gamma^{r+}_{\mathbb{AC}}=\Gamma^{ab+}_{\mathbb{CD}}=\Gamma^{r-}_{\mathbb{DB}}=\Gamma^{ab-}_{\mathbb{BA}}=\Gamma_{\mathbb{ACDBA}}\equiv\Gamma_{\circlearrowleft},\end{split}$
(21)

which implies $\Gamma_{\circlearrowright}=-\Gamma_{\circlearrowleft}$ [FIG.
2]. Using the above relations in Eq. (20), one can rewrite the steady state
energy and particle currents due to the $r$’th reservoir as

$\begin{split}J_{\rm{E}}^{r}=&\varepsilon_{\rm{R}}\Gamma_{\circlearrowleft}+(\varepsilon_{\rm{R}}+\kappa)\Gamma_{\circlearrowright}=\kappa\Gamma_{\circlearrowright};\\\
J_{\rm{N}}^{r}=&\Gamma_{\mathbb{AC}}^{r+}+\Gamma_{\mathbb{BD}}^{r+}=\Gamma_{\circlearrowleft}+\Gamma_{\circlearrowright}=0.\end{split}$
(22)

This immediately yields
$J_{\rm{Q}}^{r}=J_{\rm{E}}^{r}=\kappa\Gamma_{\circlearrowright}$. As $\rm
QD_{R}$ is only coupled with lead $r$, there will be no net particle current
driven by $r$ at the steady state, i.e., $J_{\rm{N}}^{r}=0$. However, QDL is
simultaneously coupled with leads $a$ and $b$, allowing (spin-polarized)
particle flow across $a$ and $b$ at the steady state. There is no external
source or sink associated with the system-reservoir model, so the net energy
and (spin-polarized) particle currents are conserved in the steady state. To
this end, we summarize our first set of analytical results

$\begin{split}\sum_{\lambda=a,b,r}J_{\rm{E}}^{\lambda}=0,\quad\text{and}\quad\sum_{\lambda=a,b,r}J_{\rm{N}}^{\lambda}=0,\end{split}$
(23)

which yields

$J_{\rm{E}}^{r}=-J_{\rm{E}}^{ab}=\kappa\Gamma_{\circlearrowright};\quad\quad
J_{\rm{N}}^{a}=-J_{\rm{N}}^{b}.$ (24)

The explicit form of the $\Gamma_{\circlearrowright}$ is given in Appendix-
C]. The heat current is not a conserved quantity which readily follows
$J_{\rm{Q}}^{a}+J_{\rm{Q}}^{b}\equiv J_{\rm{Q}}^{ab}\neq
J_{\rm{E}}^{ab}=-\kappa\Gamma_{\circlearrowright}$. Next, we will discuss how
the above steady-state currents are related to the entropy production rate at
near equilibrium through force-flux relation.

##  V Thermodynamic Force, Flux and Entropy Production

From the definition of the entropy change of the system [Eq. (13)], one can
write down the following relation Landi and Paternostro (2021)

$\varDelta\mathcal{S}_{\rm{s}}(t)=\Sigma(t)+\Phi(t),$ (25)

where the entropy production term $\Sigma(t)$ signifies the irreversible
contribution of the system entropy change and $\Phi(t)$ term refers to the
entropy flux, i.e., equivalent to the system entropy change due to the
reversible contribution of heat exchanges with the reservoirs. The explicit
forms of the $\Sigma(t)$ and the $\Phi(t)$ can be derived [See Appendix- D]
following Refs. Esposito  _et al._ (2010); Landi and Paternostro (2021) as:

$\begin{split}\Sigma(t)=&k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\\{\rho_{\rm{tot}}(t)\\}]\\\
&-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\left\\{\rho_{\rm{s}}(t)\prod_{\lambda}\rho_{\lambda}^{\rm{eq}}\right\\}],\\\
\Phi(t)=&k_{B}\sum_{\lambda}\Tr_{\rm{\lambda}}\left[\left\\{\rho_{\lambda}(t)-\rho_{\lambda}^{\rm{eq}}\right\\}\ln{\rho_{\lambda}^{\rm{eq}}}\right],\end{split}$
(26)

where, $\rho_{\lambda}^{\rm{eq}}$ and $\rho_{\lambda}(t)$ represent the
density operators of the $\lambda$’th reservoir under equilibrium and near
equilibrium conditions, respectively. Moreover, $\rho_{\lambda}^{\rm{eq}}$
under grand canonical ensemble, can be written as

$\rho_{\lambda}^{\rm{eq}}=\frac{e^{-\beta_{\lambda}H_{\rm{B}}^{\lambda}}}{\mathcal{Z_{\lambda}(\beta_{\lambda}})},$
(27)

where $H_{\rm{B}}^{\lambda}\equiv
H_{\rm{B}}^{\lambda\sigma}=\sum_{k}(\epsilon^{\lambda\sigma}_{k}-\mu_{\lambda\sigma})c_{{\lambda\sigma}k}^{\dagger}c_{{\lambda\sigma}k}$
and
${\mathcal{Z_{\lambda}(\beta_{\lambda}})}=\Tr[e^{-\beta_{\lambda}H_{\rm{B}}^{\lambda}}]$
is the the grand canonical partition function for the $\lambda$’th reservoir.
As a result, $\Phi(t)$ is identified as the heat exchange with the reservoir
between the final and initial time

$\Phi(t)=-k_{B}\sum_{\lambda}\beta_{\lambda}[{\langle
H_{\rm{B}}^{\lambda}\rangle}_{t}-{\langle
H_{\rm{B}}^{\lambda}\rangle}_{0}]=k_{B}\sum_{\lambda}\beta_{\lambda}\varDelta
Q_{\lambda},$ (28)

where $\varDelta Q_{\lambda}={\langle
H_{\rm{B}}^{\lambda}\rangle}_{0}-{\langle H_{\rm{B}}^{\lambda}\rangle}_{t}$
and ${\langle
H_{\rm{B}}^{\lambda}\rangle}_{t}=\Tr_{\rm{\lambda}}[\rho_{\lambda}(t)H_{\rm{B}}^{\lambda}]$.
Thus, we can rewrite Eq. (25) as

$\varDelta\mathcal{S}_{\rm{s}}(t)=\Sigma(t)+k_{B}\sum_{\lambda}\beta_{\lambda}\varDelta
Q_{\lambda},$ (29)

where we replace the second term by Eq. (28). At the steady state,
$\dot{\Sigma}$ becomes minimum Esposito  _et al._ (2010); Landi and
Paternostro (2021); Strasberg (2022) and there is no net entropy change in the
system, i.e.,
$\frac{d}{dt}\left\\{\varDelta\mathcal{S}_{\rm{s}}(t)\right\\}=0$. Therefore,
Eq. (29) yields

$\dot{\Sigma}=-k_{B}\sum_{\lambda}\beta_{\lambda}\frac{d}{dt}\left\\{\varDelta
Q_{\lambda}\right\\}=-k_{B}\sum_{\lambda}\beta_{\lambda}J_{\rm{Q}}^{\lambda}.$
(30)

Now, inserting the expression for heat current
$J^{\lambda}_{\rm{Q}}=J^{\lambda}_{\rm{E}}-\mu_{\lambda}J^{\lambda}_{\rm{N}}$
in Eq. (30), we obtain

$\dot{\Sigma}=-k_{B}\bigg{[}\beta_{a}J_{\rm{E}}^{a}+\beta_{b}J_{\rm{E}}^{b}+\beta_{r}J_{\rm{E}}^{r}\bigg{]}+k_{B}\bigg{[}\beta_{a}\mu_{a}J_{\rm{N}}^{a}+\beta_{b}\mu_{b}J_{\rm{N}}^{b}\bigg{]},$
(31)

where, we use Eq. (22) for $J_{\rm{N}}^{r}=0$. At the steady state, since, the
overall energy and particle currents are conserved, i.e.,
$J_{\rm{E}}=\sum_{\lambda}J_{\rm{E}}^{\lambda}=0$ and
$J_{\rm{N}}=\sum_{\lambda}J_{\rm{N}}^{\lambda}=0$. So, Eq. (31) reduces to

$\dot{\Sigma}=J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}+J_{\rm{E}}^{b}\mathcal{F}_{\rm{E}}^{b}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b},$
(32)

where the macroscopic thermodynamic forces conjugate to thermodynamic fluxes
$J_{\rm{E}}^{r|b}$ and $J_{\rm{N}}^{b}$ are identified as

$\mathcal{F}_{\rm{E}}^{r|b}=k_{B}(\beta_{a}-\beta_{r|b});\quad\mathcal{F}_{\rm{N}}^{b}=k_{B}(\beta_{b}\mu_{b}-\beta_{a}\mu_{a}).$
(33)

As a prerequisite, there are now two classes of reduced models that involve
coupled transport in energy and (spin-polarized) particle currents:

$\displaystyle(i)$ $\displaystyle\text{Reduced
model-I}\;(\mathcal{F}_{\rm{E}}^{b}=0):\;\;\dot{\Sigma}=J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b},$
$\displaystyle(ii)$ $\displaystyle\text{Reduced model-
II}\;(\mathcal{F}_{\rm{E}}^{r}=0):\;\;\dot{\Sigma}=J_{\rm{E}}^{b}\mathcal{F}_{\rm{E}}^{b}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b},$

where we set $\beta_{a}=\beta_{b}=\beta$ for Reduced model-I [FIG. 1b] and
$\beta_{a}=\beta_{r}=\beta^{\prime}$ for Reduced model-II [FIG. 1b]. Note that
the above two models are the only possible cases that can be derived from the
general model with pairs of thermodynamic forces and fluxes, as required for
ICC.

While, Eq. (32) links the entropy production rate to macroscopic forces and
fluxes through temperature, chemical potential, and amount of heat exchange
with the reservoirs, to grasp the underlying thermodynamic principles of
inverse current, we must identify the microscopic expressions of these
thermodynamic quantities. The ideal starting point to establish this
connection between macroscopic and microscopic thermodynamic frameworks is the
entropy production rate. Thus, rewriting the von Neumann entropy defined in
Eq. (12)

$\mathcal{S}_{\rm{s}}(t)=-k_{B}\sum_{\mathbb{i}}\rho_{\mathbb{i}}(t)\ln{\rho_{\mathbb{i}}(t)},$
(34)

in terms of the microscopic populations of the different system eigenstates,
$\\{\rho_{\mathbb{i}}\\}$, ($\mathbb{i=A,B,C,D}$), we can evaluate [See
Appendix E] the time evolution of the entropy change as Schnakenberg (1976)

$\begin{split}&\frac{d}{dt}\varDelta\mathcal{S}_{\rm{s}}(t)=k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{\rho_{\mathbb{A}}}{\rho_{\mathbb{B}}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{\rho_{\mathbb{A}}}{\rho_{\mathbb{B}}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{\rho_{\mathbb{B}}}{\rho_{\mathbb{D}}})\\\
&+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{\rho_{\mathbb{D}}}{\rho_{\mathbb{C}}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{\rho_{\mathbb{D}}}{\rho_{\mathbb{C}}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{\rho_{\mathbb{C}}}{\rho_{\mathbb{A}}})\Bigg{]},\end{split}$
(35)

where we use Eq. (7) for $\dot{\rho}_{\mathbb{i}}(t)$. By comparing the above
equation with Eq. (25), we can identify the microscopic version of entropy
production rate and entropy flux Schnakenberg (1976); Gupt  _et al._ (2024),
as follows [Appendix E]

$\begin{split}\dot{\Sigma}(t)=&k_{B}\Bigg{[}({\rm{k}}_{\mathbb{AB}}^{a+}\rho_{\mathbb{A}}-{\rm{k}}_{\mathbb{BA}}^{a-}\rho_{\mathbb{B}})\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{a-}\rho_{\mathbb{B}}})+({\rm{k}}_{\mathbb{AB}}^{b+}\rho_{\mathbb{A}}-{\rm{k}}_{\mathbb{BA}}^{b-}\rho_{\mathbb{B}})\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{b-}\rho_{\mathbb{B}}})+({\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}-{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}})\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}}{{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}}})\\\
+&({\rm{k}}_{\mathbb{DC}}^{a-}\rho_{\mathbb{D}}-{\rm{k}}_{\mathbb{CD}}^{a+}\rho_{\mathbb{C}})\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{a+}\rho_{\mathbb{C}}})+({\rm{k}}_{\mathbb{DC}}^{b-}\rho_{\mathbb{D}}-{\rm{k}}_{\mathbb{CD}}^{b+}\rho_{\mathbb{C}})\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{b+}\rho_{\mathbb{C}}})+({\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}-{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}})\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}}{{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}}})\Bigg{]},\\\
\dot{\Phi}(t)=&-k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}}{{\rm{k}}_{\mathbb{BA}}^{b-}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}}{{\rm{k}}_{\mathbb{DB}}^{r-}})+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{CD}}^{a+}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}}{{\rm{k}}_{\mathbb{CD}}^{b+}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}.\end{split}$
(36)

In Eq. (36), each term of the Schnakenberg entropy production rate
$\dot{\Sigma}$ Schnakenberg (1976), has a form of $(a-b)\ln(\frac{a}{b})$,
which guarantees the non-negativity of the entropy production rate. At the
steady state, it leads to the following form of the entropy production rate
[Appendix E]

$\begin{split}\dot{\Sigma}(t)=-\dot{\Phi}(t)=&\kappa\Gamma_{\circlearrowright}\Bigg{[}\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}+\left(\varepsilon_{\rm{L}}\Gamma^{b+}_{\mathbb{AB}}-(\varepsilon_{\rm{L}}+\kappa)\Gamma^{b-}_{\mathbb{DC}}\right)\Bigg{[}\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}})\Bigg{]}\\\
+&(\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{b-}_{\mathbb{DC}})\Bigg{[}k_{B}(1+\theta)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+k_{B}\theta\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}})\Bigg{]}\equiv
J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}+J_{\rm{E}}^{b}\mathcal{F}_{\rm{E}}^{b}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b},\end{split}$
(37)

where $\theta=(\frac{\varepsilon_{\rm{L}}}{\kappa})$ is the scaled system
parameter. The sign of $\theta$ plays a pivotal role in determining the ICC
behavior in both energy and (spin-polarized) particle currents, as will be
examined in the next section. From the above equation, all three associated
fluxes and their conjugated forces are readily identified as follows:

$\begin{split}J_{\rm{E}}^{r}=&\kappa\Gamma_{\circlearrowright}\quad;\quad
J_{\rm{N}}^{b}=(\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{b-}_{\mathbb{DC}}),\\\
J_{\rm{E}}^{b}=&\varepsilon_{\rm{L}}\Gamma_{\mathbb{AB}}^{b+}+(\varepsilon_{\rm{L}}+\kappa)\Gamma_{\mathbb{CD}}^{b+}=\varepsilon_{\rm{L}}J_{\rm{N}}^{b}+\kappa\Gamma_{\mathbb{CD}}^{b+};\end{split}$
(38)

and

$\displaystyle\mathcal{F}_{\rm{E}}^{r}=$
$\displaystyle\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}})=k_{B}(\beta_{a}-\beta_{r}),$
(39a) $\displaystyle\mathcal{F}_{\rm{E}}^{b}=$
$\displaystyle\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}})=k_{B}(\beta_{a}-\beta_{b}),$
(39b) $\displaystyle\mathcal{F}_{\rm{N}}^{b}=$ $\displaystyle
k_{B}(\theta+1)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+k_{B}\theta\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}})$
$\displaystyle=$ $\displaystyle
k_{B}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})-\varepsilon_{\rm{L}}\mathcal{F}_{\rm{E}}^{b}=k_{B}(\beta_{b}\mu_{b}-\beta_{a}\mu_{a})$
(39c)

Equation (38) is identical to the explicit expressions of the thermodynamic
fluxes obtained in Sec. IV. The equivalence between the macroscopic and
microscopic versions of the above forces [Cf. Eq. (39)] can be verified upon
using Eq. (9), followed by inserting the expressions of the FDFs. In case of
$\mathcal{F}_{\rm{E}}^{b}=0$ (i.e., $\beta_{a}=\beta_{b}$), the expression of
the particle force $\mathcal{F}_{\rm{N}}^{b}$ reduces to

$\begin{split}\mathcal{F}_{\rm{N}}^{b}=&k_{B}(\theta+1)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+k_{B}\theta\ln(\frac{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}})\\\
=&k_{B}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})=k_{B}\ln(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}})=k_{B}\beta\varDelta\mu.\end{split}$
(40)

Equations (38)-(40) are the first major results of our analysis. Equipped with
the above set of Eqs. (37)-(40), we are now ready to explore ICC in the
following section.

##  VI Result and Discussion

As outlined in Sec. V, let us analyze the reduced models with two forces and
conjugated fluxes, as a precursor of ICC. For all numerical plots, we have
considered $\gamma_{a}=\gamma_{b}=\gamma_{r}=\gamma$ and plotted appropriately
scaled dimensionless forces and fluxes. The scaled forces are taken either
zero or positive without loss of any generality.

###  VI.1 Reduced Model-I: ICC in energy current

$\S$ $\mathcal{F}_{\rm{E}}^{b}=0$ and
$\dot{\Sigma}=J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b}$:
It is intuitively clear when two forces
($\mathcal{F}_{\rm{E}}^{r},\mathcal{F}_{\rm{N}}^{b})=0$, both currents
$(J_{\rm{E}}^{r},J_{\rm{N}}^{b})$ would be zero. However, from the viewpoint
of the microscopic picture, it provides interesting insights:

If $\mathcal{F}_{\rm{N}}^{b}=0$, from Eq. (40) we can write

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}}\right)=\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)=1,$
(41)

applying Eq. (10) in the above equation, it can be shown that

$\begin{split}\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}}{{\rm{k}}_{\mathbb{BA}}^{b-}}\right)=&\left(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}}\right)=\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}}{{\rm{k}}_{\mathbb{BA}}^{ab-}}\right),\\\
\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}}{{\rm{k}}_{\mathbb{DC}}^{b-}}\right)=&\left(\frac{{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{DC}}^{a-}}\right)=\left(\frac{{\rm{k}}_{\mathbb{CD}}^{ab+}}{{\rm{k}}_{\mathbb{DC}}^{ab-}}\right).\end{split}$
(42)

Applying the above relations in Eq. (8), one obtains

$\Gamma_{\mathbb{AB}}^{ab+}=2\Gamma_{\mathbb{AB}}^{a+}=2\Gamma_{\mathbb{AB}}^{b+};\quad\quad\Gamma_{\mathbb{CD}}^{ab+}=2\Gamma_{\mathbb{CD}}^{a+}=2\Gamma_{\mathbb{CD}}^{b+},$
(43)

which when combined with Eq. (21), yields
$\Gamma_{\mathbb{AB}}^{b+}=-\Gamma_{\mathbb{CD}}^{b+}$. Following Eq. (38),
spin-polarized particle current $J_{\rm{N}}^{b}=0$, irrespective of the value
of the non-conjugate force $\mathcal{F}_{\rm{E}}^{r}$ [FIG. 3a: red line].
Thus, Eq. (33) reduces to
$\dot{\Sigma}=J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}\geq 0$. As a result,
$J_{\rm{E}}^{r}$ would always align in the direction of the
$\mathcal{F}_{\rm{E}}^{r}$, and it can be verified from the microscopic
picture as well. By considering $\beta>\beta_{r}$ we assume w.l.o.g.,
$\mathcal{F}_{\rm{E}}^{r}>0$, such that Eq. (39a) leads to the following
relation

$\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}})>0.$
(44)

Given, $\mathcal{F}_{\rm{N}}^{b}=0$, we can utilize Eq. (42) for $\kappa>0$,
to rewrite the above condition as

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{ab-}\rho_{\mathbb{B}}}\right)\left(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}}{{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}}}\right)\left(\frac{{\rm{k}}_{\mathbb{DC}}^{ab-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{ab+}\rho_{\mathbb{C}}}\right)\left(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}}{{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}}}\right)>1.$
(45)

To satisfy the above condition, at steady state, each term in the parenthesis
must be greater than $1$. Applying this in Eqs. (8) and (21), we obtain,
$\Gamma_{\circlearrowright}>0$. Similarly, for $\kappa<0$, we deduce
$\Gamma_{\circlearrowright}<0$. In both cases, the energy current is defined
as $J_{\rm{E}}^{r}=\kappa\Gamma_{\circlearrowright}>0$ [FIG. 3a: blue line].
So, ICC can occur neither in energy nor (spin-polarized) particle currents.

Figure 3: Variation of $J_{\rm{E}}^{r}$ (blue line) and $J_{\rm{N}}^{b}$ (red
line) against (a) $\mathcal{F}_{\rm{E}}^{r}$ while
$\mathcal{F}_{\rm{N}}^{b}=0$ i.e. $(\mu_{a}=\mu_{b})$ and (b)
$\mathcal{F}_{\rm{N}}^{b}$ while $\mathcal{F}_{\rm{E}}^{r}=0$ i.e.
$(\beta=\beta_{r})$. For all cases, the (dashed) dotted lines correspond to
$\kappa=(-)+1.5\hbar\gamma$, validating the condition
$|\kappa|>\varepsilon_{\rm{L}}$. Other system and bath parameters:
$\varepsilon_{\rm{L}}=1.0\hbar\gamma$, $\varepsilon_{\rm{R}}=2.5\hbar\gamma$,
$\beta_{r}=1/\hbar\gamma$ and $\mu_{b}=1.0\hbar\gamma$.

$\bullet$ When $\mathcal{F}_{\rm{E}}^{r}=0$: The non-negativity of the entropy
production rate $\dot{\Sigma}=J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b}\geq 0$,
implies that $J_{\rm{N}}^{b}$ would always align in the direction of
$\mathcal{F}_{\rm{N}}^{b}$ [FIG. 3a: red line]. This can be verified also from
the microscopic picture. Without loss of any generality, we consider
$\mathcal{F}_{\rm{N}}^{b}>0$ in Eq. (40) which yields

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}}\right)>1\quad;\quad\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)>1,$
(46)

and immediately suggests

$\begin{split}\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}}{{\rm{k}}_{\mathbb{BA}}^{b-}}\right)>&\left(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}}\right)\quad;\quad\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}}{{\rm{k}}_{\mathbb{DC}}^{b-}}\right)>\left(\frac{{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{DC}}^{a-}}\right).\end{split}$
(47)

Applying Eq. (10) in Eq. (47), it can be shown that,
${\rm{k}}_{\mathbb{BA|DC}}^{a-}>{\rm{k}}_{\mathbb{BA|DC}}^{b-}$ and
${\rm{k}}_{\mathbb{AB|CD}}^{a+}<{\rm{k}}_{\mathbb{AB|CD}}^{b+}$. Using these
relations in Eq. (8) we obtain,
$\Gamma_{\mathbb{AB|CD}}^{b+}>\Gamma_{\mathbb{AB|CD}}^{a+}$. Now, we can
define two variables $\rm{X}$ and $\rm{Y}$ as

$\rm{X}=\Gamma_{\mathbb{AB}}^{b+}-\Gamma_{\mathbb{AB}}^{a+}\quad;\quad\rm{Y}=\Gamma_{\mathbb{CD}}^{b+}-\Gamma_{\mathbb{CD}}^{a+},$
(48)

such that we can write

$\Gamma_{\mathbb{AB}}^{ab+}=2\Gamma_{\mathbb{AB}}^{b+}-\rm{X}\quad;\quad\Gamma_{\mathbb{CD}}^{ab+}=2\Gamma_{\mathbb{CD}}^{b+}-\rm{Y}.$
(49)

Using Eqs. (21) and (38), it’s straightforward to demonstrate from Eq. (49)
that under steady-state conditions,
$J_{\rm{N}}^{b}=\frac{1}{2}(\rm{X}+\rm{Y})>0$, as both $X$ and $Y$ are
positive in this case. Again, from Eq. (39a) it follows

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)=1.$
(50)

We can rearrange the above conditions as

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)=\left(\frac{{\rm{k}}_{\mathbb{BA}}^{a+}{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{DC}}^{ab-}}{{\rm{k}}_{\mathbb{AB}}^{a-}{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CD}}^{ab+}}\right),$
(51)

where the r.h.s can be re-written as

$\begin{split}\rm{r.h.s}=&\left\\{\frac{1+({{\rm{k}}_{\mathbb{AB}}^{b+}}/{{\rm{k}}_{\mathbb{AB}}^{a+}})}{1+({{\rm{k}}_{\mathbb{BA}}^{b-}/{\rm{k}}_{\mathbb{BA}}^{a-}})}\right\\}\left\\{\frac{1+({{\rm{k}}_{\mathbb{DC}}^{b-}}/{{\rm{k}}_{\mathbb{DC}}^{a-}})}{1+({{\rm{k}}_{\mathbb{CD}}^{b+}/{\rm{k}}_{\mathbb{CD}}^{a+}})}\right\\}\equiv\mathcal{P}\mathcal{Q}.\end{split}$
(52)

It is trivial to show that

$\displaystyle({{\rm{k}}_{\mathbb{AB}}^{b+}}/{{\rm{k}}_{\mathbb{AB}}^{a+}})>1\quad\text{and}\quad({{\rm{k}}_{\mathbb{BA}}^{b-}/{\rm{k}}_{\mathbb{BA}}^{a-}})<1\quad\Rightarrow\mathcal{P}>1;$
$\displaystyle({{\rm{k}}_{\mathbb{DC}}^{b-}}/{{\rm{k}}_{\mathbb{DC}}^{a-}})<1\quad\text{and}\quad({{\rm{k}}_{\mathbb{CD}}^{b+}/{\rm{k}}_{\mathbb{CD}}^{a+}})>1\quad\Rightarrow\mathcal{Q}<1,$
(53)

hence, we obtain from Eq. (52) that $\mathcal{PQ}\gtrless 1$. Further, it is
possible to verify

$\left\\{\frac{1+({{\rm{k}}_{\mathbb{AB}}^{b+}}/{{\rm{k}}_{\mathbb{AB}}^{a+}})}{1+({{\rm{k}}_{\mathbb{BA}}^{b-}/{\rm{k}}_{\mathbb{BA}}^{a-}})}\right\\}<\left\\{\frac{1+({{\rm{k}}_{\mathbb{CD}}^{b-}}/{{\rm{k}}_{\mathbb{CD}}^{a-}})}{1+({{\rm{k}}_{\mathbb{DC}}^{b+}/{\rm{k}}_{\mathbb{DC}}^{a+}})}\right\\},$
(54)

under the conditions $\mathcal{F}_{\rm{E}}^{b}\geqslant 0$ and $\kappa>0$,
which gives $\mathcal{P}\mathcal{Q}<1$. Thus, r.h.s of Eq. (51) reduces to

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)=\mathcal{P}\mathcal{Q}<1.$
(55)

For $\kappa<0$, following similar arguments, one can easily show that
$\mathcal{P}\mathcal{Q}>1$. Hence, following the same treatment done in Eqs.
(44) (45), one may conclude from Eqs. (8) and (21) that
$\Gamma_{\circlearrowright}\lessgtr 0$ for $\kappa\gtrless 0$, respectively.
This immediately suggests that the energy current,
$J_{\rm{E}}^{r}=\kappa\Gamma_{\circlearrowright}<0$, irrespective of the sign
of $\kappa$ [FIG. 3b: blue line]. Thus, we have an energy current flowing
against $\mathcal{F}_{\rm{N}}^{b}$ as pseudo-ICC, a special type of cross-
effect.

Figure 4: Variation of the energy current $J_{\rm{E}}^{r}$ with both
thermodynamic forces $\mathcal{F}_{\rm E}^{r}$ and $\mathcal{F}_{\rm N}^{b}$
for positive and negative $\kappa$: (a) $\kappa=1.5\hbar\gamma$ and (b)
$\kappa=-1.5\hbar\gamma$. For both plots, the region left to the white dashed
line signifies the genuine ICC in the energy current. Other system and bath
parameters: $\varepsilon_{\rm{L}}=1.0\hbar\gamma$,
$\varepsilon_{\rm{R}}=2.5\hbar\gamma$, $\beta_{r}=1/\hbar\gamma$, and
$\mu_{b}=1.0\hbar\gamma$.

A note on the cross-effect and the term pseudo-ICC is appropriate here. In
both cases, one force is set to zero, causing the conjugate flux to be
unaffected by its corresponding force, which lacks a specific direction. Under
these conditions, the current can flow either along or against the non-
conjugate force, as both are thermodynamically allowed. If the flux aligns
with the non-conjugate force, we call it as the normal cross-effect or simply
“cross-effect”. If the flux opposes the applied force, we classify it as a
special type of cross-effect or pseudo-ICC. While achieving genuine ICC is a
challenging task, pseudo-ICC could serve as an initial indication. In the
Reduced model-I, since pseudo-ICC behavior is observed solely in the energy
current, true ICC may only be expected in the energy current.

$\bullet$ When $(\mathcal{F}_{\rm{E}}^{r},\;\mathcal{F}_{\rm{N}}^{b})>0$:
Among the four possibilities, both fluxes could align in the direction of
their conjugate force, validating the positivity of the entropy production
rate
$\dot{\Sigma}=J_{\rm{E}}^{r}\mathcal{F}_{\rm{E}}^{r}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b}$.
However, both fluxes simultaneously can’t operate against both forces, which
would otherwise violate the laws of thermodynamics. Therefore, we can have ICC
in either $J_{\rm{E}}^{r}$ or $J_{\rm{N}}^{b}$. However, ICC in the
$J_{\rm{N}}^{b}$ is not possible in this reduced model-I, as discussed above.
Hence, we will focus on the true ICC effect in $J_{\rm{E}}^{r}$ only.

It is trivial to understand when $J_{\rm{E}}^{r}$ flows against both of the
forces, the model behaves like an autonomous refrigerator that generates
energy current against the thermal gradient. The other flux $J_{\rm{N}}^{b}$
must be positive and compensate the negative contribution of $J_{\rm{E}}^{r}$
to make the entropy production rate non-negative. Given the microscopic
condition of $\mathcal{F}_{\rm{N}}^{b}>0$, we have already verified that
$J_{\rm{N}}^{b}>0$. Again, $\mathcal{F}_{\rm E}^{r}>0$, so in view of Eqs.
(39a), the microscopic condition described by Eq. (50), is modified as

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)\gtrless
1,\quad\text{for}\quad\kappa\gtrless 0.$ (56)

Following the similar argument mentioned in Eqs. (50)- (52), we can rewrite
the above condition as

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)\gtrless\mathcal{P}\mathcal{Q}\quad\text{for}\quad\kappa\gtrless
0.$ (57)

We have already shown as long as $\mathcal{F}_{\rm{E}}^{b}\geqslant 0$, one
can write $\mathcal{P}\mathcal{Q}\lessgtr 1$ for $\kappa\gtrless 0$,
respectively. Hence, the final condition from Eq. (57) reduces to

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)\gtrless
1,$ (58)

which is not governed by the sign of $\kappa$. Now with a similar procedure
presented before, Eq. (58) can result in $\Gamma_{\circlearrowright}\gtrless
0$ for any values of $\kappa$. Hence, the energy current
$J_{\rm{E}}^{r}=\kappa\Gamma_{\circlearrowright}$, could be positive or
negative, for $\kappa\gtrless 0$. While $J_{\rm{E}}^{r}>0$ is normally
expected as both forces
$(\mathcal{F}_{\rm{E}}^{r},\;\mathcal{F}_{\rm{N}}^{b})>0$, but,
$J_{\rm{E}}^{r}<0$, signifies genuine ICC, as the energy current would flow
against both the forces, indicating the autonomous refrigeration effect, which
can be obtained in both positive and negative $\kappa$. It is evident, to
maintain a positive entropy production rate, ICC in energy current is achieved
only for a small $\mathcal{F}_{\rm E}^{r}$ and large $\mathcal{F}_{\rm
N}^{b}$, as depicted in FIG. 4. Thus, we can implement our Reduced model-I as
an autonomous refrigeration device Zhang  _et al._ (2015); Erdman  _et al._
(2018); Daré (2019); Mukherjee  _et al._ (2020); Mayrhofer  _et al._ (2021);
Barman  _et al._ (2021) i.e., a spin-thermoelectric machine that drives the
energy current from the cold reservoir to the hot reservoir by the influence
of the spin-polarized particle current. The Coefficient of performance (COP)
of such device is given by Sánchez and Büttiker (2011); Zhang  _et al._ (2023)

$\zeta=\bigg{|}\frac{J_{\rm{Q}}^{r}}{\Delta\mu\cdot J_{\rm{N}}^{b}}\bigg{|},$
(59)

where the numerator, $|J_{\rm{Q}}^{r}|=|J_{\rm{E}}^{r}|$ represents the
cooling power.

###  VI.2 Reduced Model-II: ICC in both energy and (spin-polarized) particle
currents

$\S$ $\mathcal{F}_{\rm{E}}^{r}=0$ and
$\dot{\Sigma}=J_{\rm{E}}^{b}\mathcal{F}_{\rm{E}}^{b}+J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b}$:
It is apparent that there would be no flux when both forces are zero. To
explore pseudo-ICC as a primary criterion for ICC, we proceed by setting one
non-zero force at a time. Subsequently, we will set both forces non-zero to
ascertain genuine ICC.

(i) When $\mathcal{F}_{\rm{N}}^{b}=0$ and $\mathcal{F}_{\rm{E}}^{b}>0$: Upon
setting the conditions in Eq. (39c) for the particle force, one can obtain
after a little bit of algebra

$\begin{split}&\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})=\theta\ln(\frac{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}})>0,\\\
\text{or,}&\;\left(1+\theta\right)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})=\theta\ln\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)>0,\\\
\text{or,}&\;\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})=\left(\frac{\theta}{1+\theta}\right)\ln\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)>0,\end{split}$
(60)

where,
$\frac{\theta}{1+\theta}=\frac{\varepsilon_{\rm{L}}}{\varepsilon_{\rm{L}}+\kappa}$.
So, there are two possibilities:

$\bullet$ For $\kappa>0$, it directly reduces to Eq. (46), with both $\rm{X}$
and $\rm{Y}$ positive, which finally leads to the condition
$J_{\rm{N}}^{b}>0$. Hence, the spin-induced particle current is always
positive [FIG. 5a], indicating the presence of the normal cross-effect, as the
conjugate force $\mathcal{F}_{\mathrm{N}}^{b}$ is zero.

$\bullet$ If $\kappa<0$ and $|\kappa|>\varepsilon_{\rm{L}}$, the positions of
the energy eigenstates $|\mathbb{C}\rangle$ and $|\mathbb{D}\rangle$ are
interchanged, and the Eq. (60) reduces to

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}}\right)>1\quad;\quad\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)<1.$
(61)

Following a similar procedure presented in the context of Eqs. (46) (49), one
can show that $J_{\rm{N}}^{b}=\frac{1}{2}(\rm{X}+\rm{Y})\gtrless 0$, as
$\rm{X}$ is always positive and $\rm{Y}$ is always negative following Eq.
(61). Hence, $J_{\rm{N}}^{b}$ could be either positive or negative, indicating
the possibility of both the cross-effect and pseudo-ICC [FIG. 5a], given that
the conjugate force $\mathcal{F}_{\rm{N}}^{b}=0$. As a result,
$\dot{\Sigma}=J_{\rm{E}}^{b}\mathcal{F}_{\rm{E}}^{b}$ and $J_{\rm{E}}^{b}$
aligns with the direction of $\mathcal{F}_{\rm{E}}^{b}$.

In summary, when $\mathcal{F}_{\rm{E}}^{b}$ is the only force acting on the
system, $J_{\rm{E}}^{b}$ is always positive [FIG. 5a]. If $\kappa>0$,
$J_{\rm{N}}^{b}$ is positive, indicating the cross-effect. Conversely, if
$\kappa<0$, $J_{\rm{N}}^{b}$ can be either positive or negative, suggesting
the possibility of true ICC in the spin-induced particle current. While Eq.
(46) indicates that both $|\mathbb{A}\rangle\rightarrow|\mathbb{B}\rangle$ and
$|\mathbb{C}\rangle\rightarrow|\mathbb{D}\rangle$ transitions are primarily
controlled by the same reservoir (bath $b$), Eq. (61) signifies a competition
between the two coupled reservoirs (baths $b$ and $a$) in controlling the
respective transitions. This competition between the two reservoirs results in
$J_{\rm{N}}^{b}$ being negative.

(ii) When $\mathcal{F}_{\rm{E}}^{b}=0$ and $\mathcal{F}_{\rm{N}}^{b}>0$: In
this case, $J_{\rm{N}}^{b}$ is always positive [FIG. 5b], in confirmation with
the non-negativity of the entropy production rate
$\dot{\Sigma}=J_{\rm{N}}^{b}\mathcal{F}_{\rm{N}}^{b}$. However, the energy
current $J_{\rm{E}}^{b}$ is strongly dependent on the sign of $\kappa$ and can
be expressed as $J_{\rm{E}}^{b}>\kappa\Gamma_{\mathbb{CD}}^{b+}$ [Cf.
Eq.(38)]. As demonstrated in Eqs.(50)-(55), when $\mathcal{F}_{\rm{E}}^{r}=0$,
the value of $\mathcal{P}\mathcal{Q}\lessgtr 1$ under the conditions
$\mathcal{F}_{\rm{E}}^{b}=0$ and $\kappa\gtrless 0$. This implies

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)\lessgtr
1.$ (62)

Utilizing the steady state condition [Cf. Eq. (21)], we can derive from Eq.
(62)

$\left(\frac{{\rm{k}}_{\mathbb{DC}}^{ab-}\rho_{\rm{D}}}{{\rm{k}}_{\mathbb{CD}}^{ab+}\rho_{\rm{C}}}\right)\lessgtr
1.$ (63)

The above equation suggests that $\Gamma_{\mathbb{CD}}^{ab+}\gtrless 0$.
Substituting this into Eq. (49), we obtain
$\Gamma_{\mathbb{CD}}^{b+}\gtrless\frac{\rm{Y}}{2}$ for $\kappa\gtrless 0$.
This implies that $\Gamma_{\mathbb{CD}}^{b+}$ is always positive when
$\kappa>0$, but it can be positive or negative for $\kappa<0$. Therefore, the
energy current $J_{\rm{E}}^{b}$ is always positive for $\kappa>0$, as
$J_{\rm{E}}^{b}>\kappa\Gamma_{\mathbb{CD}}^{b+}$. However, for $\kappa<0$, the
quantity $\kappa\Gamma_{\mathbb{CD}}^{b+}$ can be either positive or negative,
resulting in $J_{\rm{E}}^{b}$ being either positive or negative [FIG. 5b].

Figure 5: Variation of the currents $J_{\rm{E}}^{b}$ (bluegreen line) and
$J_{\rm{N}}^{b}$ (red line) against the thermodynamic forces (a)
$\mathcal{F}_{\rm{E}}^{b}$ while $\mathcal{F}_{\rm{N}}^{b}=0$ i.e.
$(\mu_{a}\beta^{\prime}=\mu_{b}\beta_{b})$ and (b) $\mathcal{F}_{\rm{N}}^{b}$
while $\mathcal{F}_{\rm{E}}^{b}=0$ i.e. $(\beta^{\prime}=\beta_{b})$. For all
cases, the (dashed) dotted lines correspond to $\kappa=(-)+1.5\hbar\gamma$,
validating the condition $|\kappa|>\varepsilon_{\rm{L}}$. Other system and
bath parameters: $\varepsilon_{\rm{L}}=1.0\hbar\gamma$,
$\varepsilon_{\rm{R}}=2.5\hbar\gamma$, $\beta_{b}=1/\hbar\gamma$,
$\mu_{b}=1.0\hbar\gamma$.

To summarize, when $\mathcal{F}_{\rm{N}}^{b}$ is the only non-zero force
acting on the system, the spin-induced particle current $J_{\rm{N}}^{b}$ is
always positive. The energy current $J_{\rm{E}}^{b}$ is positive only if
$\kappa>0$, while it can be either positive or negative for $\kappa<0$.
$J_{\rm{E}}^{b}>0$ represents the cross-effect, whereas $J_{\rm{E}}^{b}<0$
signifies the pseudo-ICC in the energy current, as the current flows against
the non-conjugate particle force $\mathcal{F}_{\rm{N}}^{b}$. Hence, there is a
possibility of true ICC occurring in the energy current as well within this
reduced model.

(iii) When both $(\mathcal{F}_{\rm{E}}^{b},\mathcal{F}_{\rm{N}}^{b})>0$: From
our discussion of the previous two cases, we have found that pseudo-ICC arises
in energy and particle currents only when $\kappa<0$. Thus, we will focus
solely on the $\kappa<0$ regime.

$\bullet$ ICC in particle current: Since $\mathcal{F}_{\rm{N}}^{b}>0$, it
readily follows from Eq. (39c)

$\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})>\left(\frac{{\varepsilon_{\rm{L}}}}{k_{B}}\right)\mathcal{F}_{\rm{E}}^{b}>0.$
(64)

Inserting the expression of $\mathcal{F}_{\rm{E}}^{b}$ from Eq. (39b), in the
above equation, we obtain

$\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})>\left(\frac{\varepsilon_{\rm{L}}}{\varepsilon_{\rm{L}}+\kappa}\right)\ln\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right),$
(65)

which then implies

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}}\right)>1;\quad\left(\frac{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}\right)\gtrless
1.$ (66)

So, we can continue our discussion by classifying the above relations into two
categories: (i) In the first case, both of the arguments of the above
relations are greater than 1, and (ii) in the second case, the latter argument
is considered to be less than 1. For case (i), following the same treatment
done in Eq. (46), it can be shown as a trivial case that
$J_{\rm{N}}^{b}=\frac{1}{2}(\rm{X}+\rm{Y})>0$, as both $\rm{X}$ and $\rm{Y}$
are positive, i.e., the spin-induced particle current would be positive as
long as the conjugate force $\mathcal{F}_{\rm{N}}^{b}$ is positive. For case
(ii), as the latter argument is less than 1, $\rm{Y}$ would be negative.
Again, following the same method done in case of Eq. (61), we can obtain,
$J_{\rm{N}}^{b}=\frac{1}{2}(\rm{X}+\rm{Y})\gtrless 0$, as $\rm{X}$ is positive
and $\rm{Y}$ is negative. So, we can have genuine ICC in the particle current,
defined by $J_{\rm{N}}^{b}<0$ [FIG. 6a], which indicates the flow of the spin-
induced particle current against both the non-zero positive forces
$\mathcal{F}_{\rm{N}}^{b}$ and $\mathcal{F}_{\rm{E}}^{b}$. Thus, our Reduced
model-II can be implemented as the spin-thermoelectric heat engine Donsa  _et
al._ (2014); Thierschmann  _et al._ (2015); Daré and Lombardo (2017); Singha
(2020), which drives the spin-polarized particle flux against the particle
force, influenced by the non-conjugate energy force
$\mathcal{F}_{\rm{E}}^{b}$. Finally, one can define the output power as
$|\Delta\mu\cdot J_{\rm{N}}^{b}|$, so that the efficiency of the device is
given by Sánchez and Büttiker (2011); Zhang  _et al._ (2023)

$\eta=\bigg{|}\frac{\Delta\mu\cdot J_{\rm{N}}^{b}}{J_{\rm{Q}}^{r}}\bigg{|}.$
(67)

Figure 6: Variation of (a) the spin-polarised particle current
$J_{\rm{N}}^{b}$ and (b) the energy current $J_{\rm{E}}^{r}$ with both the
thermodynamic forces for negative $\kappa$ i.e. $\kappa=-1.5\hbar\gamma$.
Genuine ICC for $J_{\rm{N}}^{b}$ occurs in the area below the white dashed
line in (a) while it occurs for $J_{\rm{E}}^{r}$ above the white dashed line
in (b). Other system and bath parameters:
$\varepsilon_{\rm{L}}=1.0\hbar\gamma$, $\varepsilon_{\rm{R}}=2.5\hbar\gamma$,
$\beta_{b}=1/\hbar\gamma$, $\mu_{b}=1.0\hbar\gamma$.

$\bullet$ ICC in energy current: Under the first case when both the arguments
of Eq. (66) are greater than 1, the spin-induced particle current
$J_{\rm{N}}^{b}>0$. As a consequence, we obtain from Eq. (38),
$J_{\rm{E}}^{b}>\kappa\Gamma_{\mathbb{CD}}^{b+}$. Again,
$\mathcal{F}_{\rm{E}}^{r}$ is anyway zero, and under the conditions
$\mathcal{F}_{\rm{E}}^{b}>0$ and $\kappa<0$, we have already shown in Eq.
(50)- (55) that $\mathcal{P}\mathcal{Q}>1$, which results in

$\left(\frac{{\rm{k}}_{\mathbb{AB}}^{ab+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{ab-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{ab-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{ab+}{\rm{k}}_{\mathbb{AC}}^{r+}}\right)>1.$
(68)

From the above relation, we have previously established in Eq. (62) that
$\Gamma_{\mathbb{CD}}^{ab+}<0$ for the $\kappa<0$, which finally leads to
$\Gamma_{\mathbb{CD}}^{b+}<\frac{\rm{Y}}{2}$, where $\rm{Y}$ is positive.
Thus, $\Gamma_{\mathbb{CD}}^{b+}$ can be both positive or negative. As a
result, energy flux $J_{\rm{E}}^{b}$ defined as
$J_{\rm{E}}^{b}>\kappa\Gamma_{\mathbb{CD}}^{b+}$, can be positive or negative,
i.e., $J_{\rm{E}}^{b}\gtrless 0$. The genuine ICC will occur when
$J_{\rm{E}}^{b}<0$ [FIG. (6)b], indicating the flow of energy current against
both non-zero positive forces $\mathcal{F}_{\rm{N}}^{b}$ and
$\mathcal{F}_{\rm{E}}^{b}$. Consequently, the Reduced model-II functions as a
thermoelectric refrigerator Zhang  _et al._ (2015); Erdman  _et al._ (2018);
Daré (2019); Mukherjee  _et al._ (2020); Mayrhofer  _et al._ (2021); Barman
_et al._ (2021), propelling energy flux against the temperature gradient,
influenced by the non-conjugate particle force $\mathcal{F}_{\rm{N}}^{b}$. The
COP for such a device is given by Eq. (59).

So, we conclude that when both forces of the Reduced model-II,
$\mathcal{F}_{\rm{E}}^{b}$ and $\mathcal{F}_{\rm{N}}^{b}$ are positive,
genuine ICC can be achieved in both energy and spin-induced particle currents.
Thus, our model can be implemented both as a spin-thermoelectric heat engine
and a refrigerator. It is evident that true ICC in the spin-induced particle
current $J_{\rm{N}}^{b}$ is obtained when both arguments of Eq. (66) are
positive and $\kappa<0$, whereas true ICC in the energy current
$J_{\rm{E}}^{b}$ is achieved, when the latter argument of Eq. (66) is
negative, under the restriction $-\kappa>\varepsilon_{\rm{L}}>0$. These
results complement the fact that two ICC regions do not overlap with each
other [FIG. 6] in compliance with the second law of thermodynamics and
positivity of the entropy production rate. It is worth mentioning that genuine
ICC can be explored in both energy and particle currents within the general
model as well, under various system-bath parameters. However, conducting
analytical studies based on microscopic descriptions will be challenging for
the general model.

##  VII Conclusions

Identification of exact thermodynamic forces and conjugate fluxes in the
presence of multiple reservoirs is crucial for any thermodynamic description
of novel phenomena. Here we present a comprehensive quantum thermodynamic
theory of the inverse current phenomenon in coupled transport, where one
induced current opposes all the thermodynamic forces present in the system.
Based on a simple variant of Sánchez-Büttiker model of a three-terminal
Coulomb-coupled quantum dots, we have examined the counter-intuitive inverse
current behavior observed in energy and spin-induced particle currents near
equilibrium situations. Our analysis focuses on the macroscopic and
microscopic correspondence of the entropy production rate, employing the
grand-canonical formalism of the Lindblad master equation and the Schnakenberg
entropy formulation. The linearity of the quantum master equation allows exact
analytical expressions for thermodynamic forces and fluxes, incorporating both
macroscopic reservoir parameters and microscopic system characteristics. It
enables us to uniquely identify all thermodynamic force-flux pairs for both
general and reduced models, facilitating a systematic analysis of genuine ICC
behavior in both energy and spin-polarized particle currents. Finally, we
illustrate that our model can function as an autonomous spin-thermoelectric
engine or refrigerator by exploiting ICC in spin-induced particle currents and
energy currents, respectively. While we demonstrate an autonomous quantum-dot
refrigerator is easier to achieve, where chemical work done by current-
carrying quantum particles assists in transferring thermal energy from a cold
to a hot bath, autonomous engine realization necessitates attractive
interaction between coupled dots. We anticipate that our findings will be
valuable in developing ICC-assisted unconventional spin-thermometric quantum
dot devices in the near future.

## Acknowledgements

AG acknowledges IITK for infrastructure and financial support. S.G.
acknowledges the Ministry of Education, Government of India, for the Prime
Minister Research Fellowship (PMRF). N.G. is thankful to CSIR for the
fellowship.

## References

  * Cleuren and den Broeck (2001) B. Cleuren and C. Van den Broeck, “Ising model for a brownian donkey,” Europhysics Letters 54, 1 (2001).
  * Cleuren and Van den Broeck (2002) B. Cleuren and C. Van den Broeck, “Random walks with absolute negative mobility,” Phys. Rev. E 65, 030101 (2002).
  * Eichhorn  _et al._ (2002) Ralf Eichhorn, Peter Reimann, and Peter Hänggi, “Brownian motion exhibiting absolute negative mobility,” Phys. Rev. Lett. 88, 190601 (2002).
  * Höpfel  _et al._ (1986) Ralph A. Höpfel, Jagdeep Shah, Peter A. Wolff, and Arthur C. Gossard, “Negative absolute mobility of minority electrons in gaas quantum wells,” Phys. Rev. Lett. 56, 2736–2739 (1986).
  * Keay  _et al._ (1995) B. J. Keay, S. Zeuner, S. J. Allen, K. D. Maranowski, A. C. Gossard, U. Bhattacharya, and M. J. W. Rodwell, “Dynamic localization, absolute negative conductance, and stimulated, multiphoton emission in sequential resonant tunneling semiconductor superlattices,” Phys. Rev. Lett. 75, 4102–4105 (1995).
  * Cleuren and Van den Broeck (2003) B. Cleuren and C. Van den Broeck, “Brownian motion with absolute negative mobility,” Phys. Rev. E 67, 055101 (2003).
  * Machura  _et al._ (2007) L. Machura, M. Kostur, P. Talkner, J. Łuczka, and P. Hänggi, “Absolute negative mobility induced by thermal equilibrium fluctuations,” Phys. Rev. Lett. 98, 040601 (2007).
  * Dandogbessi and Kenfack (2015) Bruno S. Dandogbessi and Anatole Kenfack, “Absolute negative mobility induced by potential phase modulation,” Phys. Rev. E 92, 062903 (2015).
  * Cividini  _et al._ (2018) J Cividini, D Mukamel, and H A Posch, “Driven tracer with absolute negative mobility,” Journal of Physics A: Mathematical and Theoretical 51, 085001 (2018).
  * Słapik  _et al._ (2019) A. Słapik, J. Łuczka, P. Hänggi, and J. Spiechowicz, “Tunable mass separation via negative mobility,” Phys. Rev. Lett. 122, 070602 (2019).
  * Ros  _et al._ (2005) Alexandra Ros, Ralf Eichhorn, Jan Regtmeier, Thanh Tu Duong, Peter Reimann, and Dario Anselmetti, “Absolute negative particle mobility,” Nature 436, 928–928 (2005).
  * Nagel  _et al._ (2008) J. Nagel, D. Speer, T. Gaber, A. Sterck, R. Eichhorn, P. Reimann, K. Ilin, M. Siegel, D. Koelle, and R. Kleiner, “Observation of negative absolute resistance in a josephson junction,” Phys. Rev. Lett. 100, 217001 (2008).
  * Lieb and Liniger (1963) Elliott H. Lieb and Werner Liniger, “Exact analysis of an interacting bose gas. i. the general solution and the ground state,” Phys. Rev. 130, 1605–1616 (1963).
  * Wang  _et al._ (2020) Jiao Wang, Giulio Casati, and Giuliano Benenti, “Inverse currents in hamiltonian coupled transport,” Phys. Rev. Lett. 124, 110607 (2020).
  * Zhang and Xie (2021) Yanchao Zhang and Zhenzhen Xie, “Inverse currents in coulomb-coupled quantum dots,” Phys. Rev. E 104, 064142 (2021).
  * Zhang  _et al._ (2023) Yanchao Zhang, Shuang Wang, Wanrong Li, and Mingzhuang Yue, “Inverse current induced thermoelectric conversion in a parallel-coupled double quantum dot system,” Physica Scripta 98, 105245 (2023).
  * Sánchez and Büttiker (2011) Rafael Sánchez and Markus Büttiker, “Optimal energy quanta to current conversion,” Phys. Rev. B 83, 085428 (2011).
  * Timm (2007) Carsten Timm, “Gauge theory for the rate equations: Electrodynamics on a network,” Phys. Rev. Lett. 98, 070604 (2007).
  * Walschaers  _et al._ (2013) Mattia Walschaers, Jorge Fernandez-de-Cossio Diaz, Roberto Mulet, and Andreas Buchleitner, “Optimally designed quantum transport across disordered networks,” Phys. Rev. Lett. 111, 180601 (2013).
  * Kulvelis  _et al._ (2015) Nikolaj Kulvelis, Maxim Dolgushev, and Oliver Mülken, “Universality at breakdown of quantum transport on complex networks,” Phys. Rev. Lett. 115, 120602 (2015).
  * Thingna  _et al._ (2016) Juzar Thingna, Daniel Manzano, and Jianshu Cao, “Dynamical signatures of molecular symmetries in nonequilibrium quantum transport,” Scientific Reports 6, 28027 (2016).
  * Maier  _et al._ (2019) Christine Maier, Tiff Brydges, Petar Jurcevic, Nils Trautmann, Cornelius Hempel, Ben P. Lanyon, Philipp Hauke, Rainer Blatt, and Christian F. Roos, “Environment-assisted quantum transport in a 10-qubit network,” Phys. Rev. Lett. 122, 050501 (2019).
  * Thingna  _et al._ (2020) Juzar Thingna, Daniel Manzano, and Jianshu Cao, “Magnetic field induced symmetry breaking in nonequilibrium quantum networks,” New Journal of Physics 22, 083026 (2020).
  * Ehrhardt  _et al._ (2021) Max Ehrhardt, Robert Keil, Lukas J. Maczewsky, Christoph Dittel, Matthias Heinrich, and Alexander Szameit, “Exploring complex graphs using three-dimensional quantum walks of correlated photons,” Science Advances 7, eabc5266 (2021).
  * Strasberg  _et al._ (2013) Philipp Strasberg, Gernot Schaller, Tobias Brandes, and Massimiliano Esposito, “Thermodynamics of a physical model implementing a maxwell demon,” Phys. Rev. Lett. 110, 040601 (2013).
  * Kutvonen  _et al._ (2016) Aki Kutvonen, Takahiro Sagawa, and Tapio Ala-Nissila, “Thermodynamics of information exchange between two coupled quantum dots,” Phys. Rev. E 93, 032147 (2016).
  * Ptaszyński and Esposito (2019) Krzysztof Ptaszyński and Massimiliano Esposito, “Thermodynamics of quantum information flows,” Phys. Rev. Lett. 122, 150603 (2019).
  * Esposito  _et al._ (2009) M. Esposito, K. Lindenberg, and C. Van den Broeck, “Thermoelectric efficiency at maximum power in a quantum dot,” Europhysics Letters 85, 60010 (2009).
  * Nakpathomkun  _et al._ (2010) Natthapon Nakpathomkun, H. Q. Xu, and Heiner Linke, “Thermoelectric efficiency at maximum power in low-dimensional systems,” Phys. Rev. B 82, 235428 (2010).
  * Sothmann and Büttiker (2012) Björn Sothmann and Markus Büttiker, “Magnon-driven quantum-dot heat engine,” Europhysics Letters 99, 27001 (2012).
  * Donsa  _et al._ (2014) S. Donsa, S. Andergassen, and K. Held, “Double quantum dot as a minimal thermoelectric generator,” Phys. Rev. B 89, 125103 (2014).
  * Zhang  _et al._ (2015) Yanchao Zhang, Guoxing Lin, and Jincan Chen, “Three-terminal quantum-dot refrigerators,” Phys. Rev. E 91, 052118 (2015).
  * Thierschmann  _et al._ (2015) Holger Thierschmann, Rafael Sánchez, Björn Sothmann, Fabian Arnold, Christian Heyn, Wolfgang Hansen, Hartmut Buhmann, and Laurens W. Molenkamp, “Three-terminal energy harvester with coupled quantum dots,” Nature Nanotechnology 10, 854–858 (2015).
  * Sierra  _et al._ (2016) Miguel A. Sierra, M. Saiz-Bretín, F. Domínguez-Adame, and David Sánchez, “Interactions and thermoelectric effects in a parallel-coupled double quantum dot,” Phys. Rev. B 93, 235452 (2016).
  * Whitney  _et al._ (2016) Robert S. Whitney, Rafael Sánchez, Federica Haupt, and Janine Splettstoesser, “Thermoelectricity without absorbing energy from the heat sources,” Physica E: Low-dimensional Systems and Nanostructures 75, 257–265 (2016).
  * Erdman  _et al._ (2017) Paolo Andrea Erdman, Francesco Mazza, Riccardo Bosisio, Giuliano Benenti, Rosario Fazio, and Fabio Taddei, “Thermoelectric properties of an interacting quantum dot based heat engine,” Phys. Rev. B 95, 245432 (2017).
  * Erdman  _et al._ (2018) Paolo Andrea Erdman, Bibek Bhandari, Rosario Fazio, Jukka P. Pekola, and Fabio Taddei, “Absorption refrigerators based on coulomb-coupled single-electron systems,” Phys. Rev. B 98, 045433 (2018).
  * Whitney  _et al._ (2018) Robert S. Whitney, Rafael Sánchez, and Janine Splettstoesser, “Quantum thermodynamics of nanoscale thermoelectrics and electronic devices,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by Felix Binder, Luis A. Correa, Christian Gogolin, Janet Anders, and Gerardo Adesso (Springer International Publishing, Cham, 2018) pp. 175–206.
  * Mazal  _et al._ (2019) Yair Mazal, Yigal Meir, and Yonatan Dubi, “Nonmonotonic thermoelectric currents and energy harvesting in interacting double quantum dots,” Phys. Rev. B 99, 075433 (2019).
  * Manzano  _et al._ (2020) Gonzalo Manzano, Rafael Sánchez, Ralph Silva, Géraldine Haack, Jonatan B. Brask, Nicolas Brunner, and Patrick P. Potts, “Hybrid thermal machines: Generalized thermodynamic resources for multitasking,” Phys. Rev. Res. 2, 043302 (2020).
  * Wang  _et al._ (2022) Luqin Wang, Zi Wang, Chen Wang, and Jie Ren, “Cycle flux ranking of network analysis in quantum thermal devices,” Phys. Rev. Lett. 128, 067701 (2022).
  * Ruokola and Ojanen (2011) Tomi Ruokola and Teemu Ojanen, “Single-electron heat diode: Asymmetric heat transport between electronic reservoirs through coulomb islands,” Phys. Rev. B 83, 241404 (2011).
  * Werlang  _et al._ (2014) T. Werlang, M. A. Marchiori, M. F. Cornelio, and D. Valente, “Optimal rectification in the ultrastrong coupling regime,” Phys. Rev. E 89, 062109 (2014).
  * Hartmann  _et al._ (2015) F. Hartmann, P. Pfeffer, S. Höfling, M. Kamp, and L. Worschech, “Voltage fluctuation to current converter with coulomb-coupled quantum dots,” Phys. Rev. Lett. 114, 146805 (2015).
  * Zhang  _et al._ (2017) Yanchao Zhang, Xin Zhang, Zhuolin Ye, Guoxing Lin, and Jincan Chen, “Three-terminal quantum-dot thermal management devices,” Applied Physics Letters 110, 153501 (2017).
  * Aligia  _et al._ (2020) A. A. Aligia, D. Pérez Daroca, Liliana Arrachea, and P. Roura-Bas, “Heat current across a capacitively coupled double quantum dot,” Phys. Rev. B 101, 075417 (2020).
  * Tesser  _et al._ (2022) Ludovico Tesser, Bibek Bhandari, Paolo Andrea Erdman, Elisabetta Paladino, Rosario Fazio, and Fabio Taddei, “Heat rectification through single and coupled quantum dots,” New Journal of Physics 24, 035001 (2022).
  * Ghosh  _et al._ (2022) Shuvadip Ghosh, Nikhil Gupt, and Arnab Ghosh, “Universal behavior of the coulomb-coupled fermionic thermal diode,” Entropy 24 (2022).
  * Esposito  _et al._ (2012) Massimiliano Esposito, Niraj Kumar, Katja Lindenberg, and Christian Van den Broeck, “Stochastically driven single-level quantum dot: A nanoscale finite-time thermodynamic machine and its various operational modes,” Phys. Rev. E 85, 031117 (2012).
  * Jiang  _et al._ (2015) Jian-Hua Jiang, Manas Kulkarni, Dvira Segal, and Yoseph Imry, “Phonon thermoelectric transistors and rectifiers,” Phys. Rev. B 92, 045309 (2015).
  * Yang  _et al._ (2019) Jing Yang, Cyril Elouard, Janine Splettstoesser, Björn Sothmann, Rafael Sánchez, and Andrew N. Jordan, “Thermal transistor and thermometer based on coulomb-coupled conductors,” Phys. Rev. B 100, 045418 (2019).
  * Ghosh  _et al._ (2012) Arnab Ghosh, Sudarson Sekhar Sinha, and Deb Shankar Ray, “Fermionic oscillator in a fermionic bath,” Phys. Rev. E 86, 011138 (2012).
  * Zhang  _et al._ (2018) Yanchao Zhang, Zhimin Yang, Xin Zhang, Bihong Lin, Guoxing Lin, and Jincan Chen, “Coulomb-coupled quantum-dot thermal transistors,” EPL (Europhysics Letters) 122, 17002 (2018).
  * Agrawal and Solomon (2004) Mukul Agrawal and G. S. Solomon, “Quantum-well band structure effects on the emission polarization from a spin-polarized electron reservoir,” Applied Physics Letters 85, 1820–1822 (2004), https://pubs.aip.org/aip/apl/article-pdf/85/10/1820/19329056/1820_1_online.pdf .
  * Vandaele  _et al._ (2017) Koen Vandaele, Sarah J. Watzman, Benedetta Flebus, Arati Prakash, Yuanhua Zheng, Stephen R. Boona, and Joseph P. Heremans, “Thermal spin transport and energy conversion,” Materials Today Physics 1, 39–49 (2017).
  * Gupt  _et al._ (2021) Nikhil Gupt, Srijan Bhattacharyya, and Arnab Ghosh, “Statistical generalization of regenerative bosonic and fermionic stirling cycles,” Phys. Rev. E 104, 054130 (2021).
  * Gupt  _et al._ (2023) Nikhil Gupt, Shuvadip Ghosh, and Arnab Ghosh, “Top-ranked cycle flux network analysis of molecular photocells,” Phys. Rev. E 108, 034305 (2023).
  * Gupt  _et al._ (2024) Nikhil Gupt, Shuvadip Ghosh, and Arnab Ghosh, “Graph theoretic analysis of three-terminal quantum dot thermocouples: Onsager relations and spin-thermoelectric effects,” Phys. Rev. B 109, 125124 (2024).
  * Breuer and Petruccione (2007) Heinz-Peter Breuer and Francesco Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Strasberg (2022) Philipp Strasberg,  _Quantum Stochastic Thermodynamics: Foundations and Selected Applications_ (Oxford University Press, 2022).
  * Schnakenberg (1976) J. Schnakenberg, “Network theory of microscopic and macroscopic behavior of master equation systems,” Rev. Mod. Phys. 48, 571–585 (1976).
  * Esposito  _et al._ (2010) Massimiliano Esposito, Katja Lindenberg, and Christian Van den Broeck, “Entropy production as correlation between system and reservoir,” New Journal of Physics 12, 013013 (2010).
  * Landi and Paternostro (2021) Gabriel T. Landi and Mauro Paternostro, “Irreversible entropy production: From classical to quantum,” Rev. Mod. Phys. 93, 035008 (2021).
  * Kosloff (2013) Ronnie Kosloff, “Quantum thermodynamics: A dynamical viewpoint,” Entropy 15, 2100–2128 (2013).
  * Daré (2019) A.-M. Daré, “Comparative study of heat-driven and power-driven refrigerators with coulomb-coupled quantum dots,” Phys. Rev. B 100, 195427 (2019).
  * Mukherjee  _et al._ (2020) Swarnadip Mukherjee, Bitan De, and Bhaskaran Muralidharan, “Three-terminal vibron-coupled hybrid quantum dot thermoelectric refrigeration,” Journal of Applied Physics 128, 234303 (2020).
  * Mayrhofer  _et al._ (2021) R. David Mayrhofer, Cyril Elouard, Janine Splettstoesser, and Andrew N. Jordan, “Stochastic thermodynamic cycles of a mesoscopic thermoelectric engine,” Phys. Rev. B 103, 075404 (2021).
  * Barman  _et al._ (2021) Anamika Barman, Surojit Halder, Shailendra K. Varshney, Gourab Dutta, and Aniket Singha, “Realistic nonlocal refrigeration engine based on coulomb-coupled systems,” Phys. Rev. E 103, 012131 (2021).
  * Daré and Lombardo (2017) A.-M. Daré and P. Lombardo, “Powerful coulomb-drag thermoelectric engine,” Phys. Rev. B 96, 115414 (2017).
  * Singha (2020) Aniket Singha, “A realistic non-local heat engine based on Coulomb-coupled systems,” Journal of Applied Physics 127, 234903 (2020).
  * Gupt  _et al._ (2022) Nikhil Gupt, Srijan Bhattacharyya, Bikash Das, Subhadeep Datta, Victor Mukherjee, and Arnab Ghosh, “Floquet quantum thermal transistor,” Phys. Rev. E 106, 024110 (2022).

##  Appendix A Derivation of the LME

We present the interaction picture master equation to study the system’s
dynamical evolution that exchanges particles and energy with the reservoirs.
Let us start with the tunneling Hamiltonian

$\begin{split}H_{\rm{T}}=H_{\rm{T}}^{\rm{L}\downarrow
a}+H_{\rm{T}}^{\rm{L}\downarrow b}+H_{\rm{T}}^{\rm{R}\uparrow r}\end{split}$
(69)

where, $H_{\rm{T}}^{\rm{L}\downarrow a(b)}=\hbar\sum_{k}[t^{\rm{L}\downarrow
a(b)}_{k}c_{a(b)\downarrow
k}^{\dagger}d_{{\rm{L}}\downarrow}+t^{\rm{L}\downarrow
a(b)*}_{k}d_{{\rm{L}}\downarrow}^{\dagger}c_{a(b)\downarrow k}]$ and
$H_{\rm{T}}^{\rm{R}\uparrow r}=\hbar\sum_{k}[t^{\rm{R}\uparrow
r}_{k}c_{r\uparrow k}^{\dagger}d_{{\rm{R}}\uparrow}+t^{\rm{R}\uparrow
r*}_{k}d_{{\rm{R}}\uparrow}^{\dagger}c_{r\uparrow k}]$. To derive the master
equation, we start with the von Neumann equation for the total density matrix
$\rho_{\rm{tot}}$ in the interaction picture

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

Integrating the above equation, and taking a trace over the bath degrees of
freedom, one obtains

$\begin{split}\frac{\partial}{\partial
t}\rho_{\rm{s}}(t)=\frac{1}{(i\hbar)^{2}}\int_{0}^{t}dt^{\prime}{\rm
Tr}_{a,b,r}[H_{\rm{T}}(t),[H_{\rm{T}}(t-t^{\prime}),\rho_{\rm{tot}}(t^{\prime})]],\end{split}$
(71)

where, ${\rm Tr}_{a,b,r}$ refers to the trace over each bath degrees of
freedom and ${\rm Tr}_{a,b,r}\\{\rho_{\rm{tot}}(t)\\}=\rho_{\rm{s}}(t)$
denotes the reduced density operator for the system. We also assume that ${\rm
Tr}_{a,b,r}[H_{\rm{T}}(t),\rho_{\rm{tot}}(0)]=0$. Under the Born-Markov
approximation, the above equation can be rewritten as Breuer and Petruccione
(2007); Gupt  _et al._ (2022); Ghosh  _et al._ (2022)

$\begin{split}\dot{\rho_{\rm{s}}}(t)=\frac{1}{(i\hbar)^{2}}\sum_{\lambda=a,b,r}\int_{0}^{\infty}dt^{\prime}\
{\rm
Tr}_{a,b,r}[H_{\rm{T}}^{\lambda}(t),[H_{\rm{T}}^{\lambda}(t-t^{\prime}),\rho_{\rm{s}}(t)\otimes\rho_{a}\otimes\rho_{b}\otimes\rho_{r}]],\end{split}$
(72)

where we use the following properties of the bath operators ${\rm
Tr}_{\lambda}\\{{c}_{\lambda}(t)\rho_{\lambda}\\}=0={\rm
Tr}_{\lambda}\\{{c}^{\dagger}_{\lambda}(t)\rho_{\lambda}\\}$ and ${\rm
Tr}_{a,b,r}\\{[H^{\lambda}_{\rm{T}}(t),[H^{\nu}_{\rm{T}}(t-t^{\prime}),{\rho_{\rm{s}}}(t)\otimes\rho_{a}\otimes\rho_{b}\otimes\rho_{r}]]\\}=0;\lambda\neq\nu;\lambda,\nu={a},{b},{r}$.
Now, in the above equation, we use the interaction picture system and bath
operators

$\displaystyle d_{\alpha\sigma}(t)$ $\displaystyle=$ $\displaystyle
e^{{iH_{\rm{s}}t}/{\hbar}}d_{\alpha\sigma}e^{{-iH_{\rm{s}}t}/{\hbar}}=\sum_{\omega_{\mathbb{ij}}>0}e^{{{-i{\omega_{\mathbb{ij}}}}t}/{\hbar}}d_{\alpha\sigma};\quad\quad{\alpha\sigma=\rm{L}\downarrow,\rm{R}\uparrow}$
$\displaystyle c_{\lambda\sigma}(t)$ $\displaystyle=$ $\displaystyle
e^{{iH_{\rm{B}}t}/{\hbar}}c_{\lambda\sigma}e^{{-iH_{\rm{B}}t}/{\hbar}}=\sum_{k}e^{{{-i({\epsilon_{k}^{\lambda\sigma}-\mu_{\lambda}})t}/{\hbar}}}c_{\lambda\sigma};\quad\quad{\lambda\sigma=a\downarrow,b\downarrow,r\uparrow}$
(73)

and their hermitian adjoints, where $\omega_{\mathbb{ij}}$ is defined as the
transition energy for the transition between the system eigenstates
$|\mathbb{i}\rangle$ and $|\mathbb{j}\rangle$. Eliminating the high-frequency
oscillating terms by the standard procedure of secular approximation, one can
finally derive the master equation in the following form

$\dot{\rho}_{\rm{s}}(t)=\mathcal{L}_{a}[\rho_{\rm{s}}(t)]+\mathcal{L}_{b}[\rho_{\rm{s}}(t)]+\mathcal{L}_{r}[\rho_{\rm{s}}(t)],$
(74)

where the Lindblad operators $\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]$ are
given by

$\displaystyle\mathcal{L}_{\lambda}[\rho_{\rm{s}}(t)]=\sum_{\\{\omega_{\alpha}\\}>0}$
$\displaystyle\mathcal{G}_{{\lambda}}(\omega_{\alpha})$
$\displaystyle\left[d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})\rho_{\rm{s}}d_{\alpha\sigma}(\omega_{\alpha})-\frac{1}{2}{\\{\rho_{\rm{s}}},d_{\alpha\sigma}(\omega_{\alpha})d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})\\}\right]$
$\displaystyle+$ $\displaystyle\mathcal{G}_{{\lambda}}(-\omega_{\alpha})$
$\displaystyle\left[d_{\alpha\sigma}(\omega_{\alpha})\rho_{\rm{s}}d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})-\frac{1}{2}{\\{\rho_{\rm{s}}},d_{\alpha\sigma}^{\dagger}(\omega_{\alpha})d_{\alpha\sigma}(\omega_{\alpha})\\}\right].$

In the above equation, we define the temperature-dependent bath spectral
functions as

$\mathcal{G}_{{\lambda}}(\omega_{\alpha})=\gamma_{\lambda}(\omega_{\alpha})f^{+}_{\lambda}(\omega_{\alpha});\quad\mathcal{G}_{{\lambda}}(-\omega_{\alpha})=\gamma_{\lambda}(\omega_{\alpha})f^{-}_{\lambda}(\omega_{\alpha}).$
(76)

The function $f^{\pm}_{\lambda}(\omega_{\mathbb{ij}})$ represents the Fermi
distribution functions (FDF) which are obtained by tracing over the bath
density operator, for example,
$f_{\lambda\sigma}^{+}\big{(}\omega_{\mathbb{ij}}\big{)}=\Tr_{\lambda\sigma}\big{(}c^{\dagger}_{\lambda\sigma}c_{\lambda\sigma}\rho_{\lambda\sigma}\big{)}$,
and
$f_{\lambda\sigma}^{-}\big{(}\omega_{\mathbb{ij}}\big{)}=\Tr_{\lambda\sigma}\big{(}c_{\lambda\sigma}c^{\dagger}_{\lambda\sigma}\rho_{\lambda\sigma}\big{)}$,
where the bath operators $c^{\dagger}_{\lambda\sigma}$ and $c_{\lambda\sigma}$
obey anti-commutation relation and the reservoir governs the transition
between eigenstate $|\mathbb{i}\rangle$ to $|\mathbb{j}\rangle$, that costs
$\omega_{\mathbb{ij}}$ amount of energy. The explicit expressions of the FDFs
for various transitions are listed below:

$\begin{split}f^{+}_{a|b}(\omega_{\mathbb{AB}})=\left[{1+\exp\bigg({\frac{\varepsilon_{\rm{L}}-\mu_{a|b}}{k_{B}T_{a|b}}}\bigg{missing})}\right]^{-1}\equiv&f_{a|b}^{1+}\quad,\quad
f^{+}_{a|b}(\omega_{\mathbb{CD}})=\left[{1+\exp\bigg({\frac{\varepsilon_{\rm{L}}+\kappa-\mu_{a|b}}{k_{B}T_{a|b}}}\bigg{missing})}\right]^{-1}\equiv
f_{a|b}^{2+},\\\
f^{-}_{a|b}(\omega_{\mathbb{BA}})=\left[{1+\exp\bigg({\frac{-\varepsilon_{\rm{L}}+\mu_{a|b}}{k_{B}T_{a|b}}}\bigg{missing})}\right]^{-1}\equiv&f_{a|b}^{1-}\quad,\quad
f^{-}_{a|b}(\omega_{\mathbb{CD}})=\left[{1+\exp\bigg({\frac{-\varepsilon_{\rm{L}}-\kappa+\mu_{a|b}}{k_{B}T_{a|b}}}\bigg{missing})}\right]^{-1}\equiv
f_{a|b}^{2-},\\\
f^{+}_{r}(\omega_{\mathbb{AC}})=\left[{1+\exp\bigg({\frac{\varepsilon_{\rm{R}}-\mu_{r}}{k_{B}T_{r}}}\bigg{missing})}\right]^{-1}\equiv&f_{r}^{1+}\quad,\quad
f^{+}_{r}(\omega_{\mathbb{BD}})=\left[{1+\exp\bigg({\frac{\varepsilon_{\rm{R}}+\kappa-\mu_{r}}{k_{B}T_{r}}}\bigg{missing})}\right]^{-1}\equiv
f_{r}^{2+},\\\
f^{-}_{r}(\omega_{\mathbb{CA}})=\left[{1+\exp\bigg({\frac{-\varepsilon_{\rm{R}}+\mu_{r}}{k_{B}T_{r}}}\bigg{missing})}\right]^{-1}\equiv&f_{r}^{1-}\quad,\quad
f^{-}_{r}(\omega_{\mathbb{DB}})=\left[{1+\exp\bigg({\frac{-\varepsilon_{\rm{R}}-\kappa+\mu_{r}}{k_{B}T_{r}}}\bigg{missing})}\right]^{-1}\equiv
f_{r}^{2-}.\end{split}$ (77)

##  Appendix B Explicit expressions of various currents

From Eq. (20), one can evaluate the explicit expressions of both energy and
particle currents for each reservoir in the following way

$\begin{split}J^{a}_{\rm{E}}=&\omega_{\mathbb{AB}}\Gamma_{\mathbb{AB}}^{a+}+\omega_{\mathbb{CD}}\Gamma_{\mathbb{CD}}^{a+}=\varepsilon_{\rm{L}}\Gamma_{\mathbb{AB}}^{a+}+(\varepsilon_{\rm{L}}+\kappa)\Gamma_{\mathbb{CD}}^{a+};\quad
J_{\rm{E}}^{b}=\omega_{\mathbb{AB}}\Gamma_{\mathbb{AB}}^{b+}+\omega_{\mathbb{CD}}\Gamma_{\mathbb{CD}}^{b+}=\varepsilon_{\rm{L}}\Gamma_{\mathbb{AB}}^{b+}+(\varepsilon_{\rm{L}}+\kappa)\Gamma_{\mathbb{CD}}^{b+};\\\
J_{\rm{E}}^{r}=&\omega_{\mathbb{AC}}\Gamma_{\mathbb{AC}}^{r+}+\omega_{\mathbb{BD}}\Gamma_{\mathbb{BD}}^{r+}=\varepsilon_{\rm{R}}\Gamma_{\mathbb{AC}}^{r+}+(\varepsilon_{\rm{R}}+\kappa)\Gamma_{\mathbb{BD}}^{r+};\quad
J_{\rm{N}}^{a}=\Gamma_{\mathbb{AB}}^{a+}+\Gamma_{\mathbb{CD}}^{a+};J_{\rm{N}}^{b}=\Gamma_{\mathbb{AB}}^{b+}+\Gamma_{\mathbb{CD}}^{b+};J_{\rm{N}}^{r}=\Gamma_{\mathbb{AC}}^{r+}+\Gamma_{\mathbb{BD}}^{r+}.\end{split}$
(78)

Inserting the above relations, the expression of the heat current associated
with each reservoir can be evaluated as

$\begin{split}J_{\rm{Q}}^{a}=&(\varepsilon_{\rm{L}}-\mu_{a})\Gamma_{\mathbb{AB}}^{a+}+(\varepsilon_{\rm{L}}+\kappa-\mu_{a})\Gamma_{\mathbb{CD}}^{a+};\quad
J_{\rm{Q}}^{b}=(\varepsilon_{\rm{L}}-\mu_{b})\Gamma_{\mathbb{AB}}^{b+}+(\varepsilon_{\rm{L}}+\kappa-\mu_{b})\Gamma_{\mathbb{CD}}^{b+};\quad
J_{\rm{Q}}^{r}=(\varepsilon_{\rm{R}}-\mu_{r})\Gamma_{\mathbb{AC}}^{r+}+(\varepsilon_{\rm{R}}+\kappa-\mu_{r})\Gamma_{\mathbb{BD}}^{r+}.\end{split}$
(79)

##  Appendix C Expression of the steady state transition rate

To determine the steady-state transition rate, we rewrite Eq. (7)-(9) as

$\mathcal{M}\begin{bmatrix}\rho_{\mathbb{A}}\\\ \rho_{\mathbb{B}}\\\
\rho_{\mathbb{C}}\\\ \rho_{\mathbb{D}}\end{bmatrix}=\begin{bmatrix}0\\\ 0\\\
0\\\ 1\end{bmatrix},$ (80)

subject to the condition
$\rho_{\mathbb{A}}+\rho_{\mathbb{B}}+\rho_{\mathbb{C}}+\rho_{\mathbb{D}}=1$,
and

$\begin{split}\mathcal{M}=\begin{bmatrix}-f_{r}^{1+}-f_{a}^{1+}-f_{b}^{1+}&{f}_{a}^{1-}+{f}_{b}^{1-}&f_{r}^{1+}&0\\\
f_{a}^{1+}+f_{b}^{1+}&-{f}_{a}^{1-}-{f}_{b}^{1-}-f_{r}^{2+}&0&{f}_{r}^{2-}\\\
f_{r}^{1+}&0&-f_{a}^{2+}-f_{b}^{2+}-{f}_{r}^{1-}&{f}_{a}^{2-}+{f}_{b}^{2-}\\\
1&1&1&1\end{bmatrix},\end{split}$ (81)

where, for the sake of simplicity of our analysis, we assume that
$\gamma_{a}\simeq\gamma_{b}\simeq\gamma_{r}\equiv\gamma$. Solving Eq. (80)
with the above matrix $\mathcal{M}$, one can obtain the steady state
population $\\{\rho_{\mathbb{i}}\\}$ in terms of which we can evaluate the
explicit expression

$\Gamma_{\circlearrowright}=-\Gamma_{\circlearrowleft}=\gamma\Bigg{[}\frac{f_{ab}^{1+}[f_{ab}^{2+}(f_{r}^{2+}-f_{r}^{1+})+2f_{r}^{2+}(f_{r}^{1+}-1)]-2f_{r}^{1+}f_{ab}^{2+}(f_{r}^{2+}-1)}{3f_{ab}^{1+}(f_{r}^{1+}-f_{r}^{2+})-6+3f_{ab}^{2+}(f_{r}^{2+}-f_{r}^{1+})}\Bigg{]},$
(82)

where, we define $f_{ab}^{1+(2+)}=f_{a}^{1+(2+)}+f_{b}^{1+(2+)}$.

##  Appendix D Non-negativity of entropy production rate

One can evaluate the expression of the entropy production ($\Sigma$) from the
entropy change of the system
($\varDelta\mathcal{S}_{\rm{s}}=\mathcal{S}_{\rm{s}}(t)-\mathcal{S}_{\rm{s}}(0)$),
which is defined as

$\begin{split}\varDelta\mathcal{S}_{\rm{s}}(t)=-k_{B}\Tr_{\rm{s}}[\rho_{\rm{s}}(t)\ln{\rho_{\rm{s}}(t)}]+k_{B}\Tr_{\rm{s}}[\rho_{\rm{s}}(0)\ln{\rho_{\rm{s}}(0)}]=-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln{\rho_{\rm{s}}(t)}]+k_{B}\Tr[\rho_{\rm{tot}}(0)\ln{\rho_{\rm{s}}(0)}].\end{split}$
(83)

We assume that the initial equilibrium state, $\rho_{\rm{tot}}(0)$ does not
display any entanglement or correlation between the system and the
environment. Therefore

$\rho_{\rm{tot}}(0)=\rho_{\rm{s}}(0)\prod_{\lambda}\rho_{\rm{\lambda}}^{\rm{eq}}.$
(84)

Inserting Eq. (84) into Eq. (83), we can continue as follows

$\begin{split}\varDelta\mathcal{S}_{\rm{s}}(t)=&-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln{\rho_{\rm{s}}(t)}]+k_{B}\Tr[\rho_{\rm{tot}}(0)\ln{\rho_{\rm{tot}}(0)}]-k_{B}\sum_{\lambda}\Tr[\rho_{\rm{tot}}(0)\ln{\rho_{\lambda}^{\rm{eq}}}]\\\
=&-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\left\\{{\rho_{\rm{s}}(t)}\prod_{\lambda}\rho_{\lambda}^{\rm{eq}}\right\\}]+k_{B}\Tr[\rho_{\rm{tot}}(0)\ln{\rho_{\rm{tot}}(0)}]+k_{B}\sum_{\lambda}\Tr_{\lambda}[\\{\rho_{\lambda}(t)-\rho_{\lambda}^{\rm{eq}}\\}\ln{\rho_{\lambda}^{\rm{eq}}}].\end{split}$
(85)

Again, $\rho_{\rm{tot}}(t)$ and $\rho_{\rm{tot}}(0)$ are related through
unitary evolution
$\rho_{\rm{tot}}(t)=\mathcal{U}\rho_{\rm{tot}}(0)\mathcal{U}^{\dagger}$, which
implies
$\Tr[\rho_{\rm{tot}}(t)\ln\rho_{\rm{tot}}(t)]=\Tr[\rho_{\rm{tot}}(0)\ln{\rho_{\rm{tot}}(0)}]$.
Applying this relation in the above equation, the final expression of the
entropy change of the system can be calculated as

$\varDelta\mathcal{S}_{\rm{s}}(t)=-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\left\\{{\rho_{\rm{s}}(t)}\prod_{\lambda}\rho_{\lambda}^{\rm{eq}}\right\\}]+k_{B}\Tr[\rho_{\rm{tot}}(t)\ln{\rho_{\rm{tot}}(t)}]+k_{B}\sum_{\lambda}\Tr_{\lambda}[\\{\rho_{\lambda}(t)-\rho_{\lambda}^{\rm{eq}}\\}\ln{\rho_{\lambda}^{\rm{eq}}}].$
(86)

The last term of the above equation can be identified as the entropy flow
($\Phi$), representing the reversible contribution to the system entropy
change due to heat exchange with the reservoirs. A comparison of the above
equation with Eq. (25), defines the entropy production, representing the
irreversible contribution to the entropy change of the system and the entropy
flow as

$\Sigma(t)=k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\\{\rho_{\rm{tot}}(t)\\}]-k_{B}\Tr[\rho_{\rm{tot}}(t)\ln\left\\{\rho_{\rm{s}}(t)\prod_{\lambda}\rho_{\lambda}^{\rm{eq}}\right\\}];\qquad\Phi(t)=k_{B}\sum_{\lambda}\Tr_{\rm{\lambda}}\left[\left\\{\rho_{\lambda}(t)-\rho_{\lambda}^{\rm{eq}}\right\\}\ln{\rho_{\lambda}^{\rm{eq}}}\right].$
(87)

The above expression of the entropy production can be expressed in terms of
the relative entropy

$\Sigma(t)\equiv\mathcal{D}\Bigg{[}\rho_{\rm{tot}}(t)||\left\\{\rho_{\rm{s}}(t)\prod_{\lambda}\rho_{\lambda}^{\rm{eq}}\right\\}\Bigg{]},$
(88)

where $\mathcal{D}[\rho||{\rho^{\prime}}]$ is the quantum relative entropy
between two density matrices $\rho$ and ${\rho^{\prime}}$, defined via

$\mathcal{D}[\rho||{\rho^{\prime}}]:=\Tr[\rho\ln\rho]-\Tr[\rho\ln{\rho^{\prime}}].$
(89)

The non-negativity of relative entropy is affirmed, attaining a value of zero
solely in the case of complete matrix identity. While the non-negativity of
entropy production doesn’t imply the same for its rate, in the limit of large
reservoirs, $\Sigma(t)$ is expected to converge to a convex, monotonically
increasing function of time Esposito  _et al._ (2010). In the same limit, if
the system dynamics are described by a Markovian quantum Lindblad master
equation, implying entropy production as a convex functional of the system
density matrix Breuer and Petruccione (2007); Esposito  _et al._ (2010), the
entropy production rate $\dot{\Sigma}(t)$ would eventually be positive, only
reaching zero for the equilibrium state.

##  Appendix E Derivation of the microscopic definition of the entropy
production rate

Considering the von Neumann entropy defined as
$\mathcal{S}_{\rm{s}}(t)=-k_{B}\sum_{\mathbb{i}}\rho_{\mathbb{i}}(t)\ln{\rho_{\mathbb{i}}(t)}$,
where, $\rho_{\mathbb{i}}$ signifies populations of the system eigenstates
($\mathbb{i=A,B,C,D}$), the change in the system entropy is given by

$\begin{split}\varDelta\mathcal{S}_{\rm{s}}(t)=&\mathcal{S}_{\rm{s}}(t)-\mathcal{S}_{\rm{s}}(0)=-k_{B}\sum_{\mathbb{i}}\rho_{\mathbb{i}}(t)\ln{\rho_{\mathbb{i}}(t)}+k_{B}\sum_{\mathbb{i}}\rho_{\mathbb{i}}(0)\ln{\rho_{\mathbb{i}}(0)}.\end{split}$
(90)

So, the time evolution of the entropy change can be evaluated as

$\begin{split}\frac{d}{dt}\varDelta\mathcal{S}_{\rm{s}}(t)=&-k_{B}\sum_{\mathbb{i}}\dot{\rho_{\mathbb{i}}}(t)\ln{\rho_{\mathbb{i}}(t)}+k_{B}\sum_{\mathbb{i}}\dot{\rho_{\mathbb{i}}}(0)\ln{\rho_{\mathbb{i}}(0)}=-k_{B}\sum_{\mathbb{i}}\dot{\rho_{\mathbb{i}}}(t)\ln{\rho_{\mathbb{i}}(t)}\equiv-
k_{B}\sum_{\mathbb{i}}\dot{\rho_{\mathbb{i}}}\ln{\rho_{\mathbb{i}}}.\end{split}$
(91)

Now, using Eq. (7) for $\dot{\rho_{\mathbb{i}}}$, we recover Eq. (35) of the
main text:

$\begin{split}\frac{d}{dt}\varDelta\mathcal{S}_{\rm{s}}(t)=k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{\rho_{\mathbb{A}}}{\rho_{\mathbb{B}}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{\rho_{\mathbb{A}}}{\rho_{\mathbb{B}}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{\rho_{\mathbb{B}}}{\rho_{\mathbb{D}}})+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{\rho_{\mathbb{D}}}{\rho_{\mathbb{C}}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{\rho_{\mathbb{D}}}{\rho_{\mathbb{C}}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{\rho_{\mathbb{C}}}{\rho_{\mathbb{A}}})\Bigg{]}.\end{split}$
(92)

By comparing the above equation with Eq. (25), we obtain Eq. (36) for the
general expressions of $\dot{\Sigma}(t)$ and $\dot{\Phi}(t)$:

$\begin{split}\dot{\Sigma}(t)=&k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{a-}\rho_{\mathbb{B}}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{b-}\rho_{\mathbb{B}}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}}{{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}}})+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{a+}\rho_{\mathbb{C}}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{b+}\rho_{\mathbb{C}}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}}{{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}}})\Bigg{]}\\\
=&k_{B}\Bigg{[}({\rm{k}}_{\mathbb{AB}}^{a+}\rho_{\mathbb{A}}-{\rm{k}}_{\mathbb{BA}}^{a-}\rho_{\mathbb{B}})\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{a-}\rho_{\mathbb{B}}})+({\rm{k}}_{\mathbb{AB}}^{b+}\rho_{\mathbb{A}}-{\rm{k}}_{\mathbb{BA}}^{b-}\rho_{\mathbb{B}})\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}\rho_{\mathbb{A}}}{{\rm{k}}_{\mathbb{BA}}^{b-}\rho_{\mathbb{B}}})+({\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}-{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}})\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}\rho_{\mathbb{B}}}{{\rm{k}}_{\mathbb{DB}}^{r-}\rho_{\mathbb{D}}})\\\
+&({\rm{k}}_{\mathbb{DC}}^{a-}\rho_{\mathbb{D}}-{\rm{k}}_{\mathbb{CD}}^{a+}\rho_{\mathbb{C}})\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{a+}\rho_{\mathbb{C}}})+({\rm{k}}_{\mathbb{DC}}^{b-}\rho_{\mathbb{D}}-{\rm{k}}_{\mathbb{CD}}^{b+}\rho_{\mathbb{C}})\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}\rho_{\mathbb{D}}}{{\rm{k}}_{\mathbb{CD}}^{b+}\rho_{\mathbb{C}}})+({\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}-{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}})\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}\rho_{\mathbb{C}}}{{\rm{k}}_{\mathbb{AC}}^{r+}\rho_{\mathbb{A}}})\Bigg{]},\\\
\dot{\Phi}(t)=&-k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}}{{\rm{k}}_{\mathbb{BA}}^{b-}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}}{{\rm{k}}_{\mathbb{DB}}^{r-}})+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{CD}}^{a+}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}}{{\rm{k}}_{\mathbb{CD}}^{b+}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}.\end{split}$
(93)

There is no net entropy change in the system at the steady state, which
reduces Eq. (37) for the form of the entropy production rate:

$\begin{split}\dot{\Sigma}(t)&=-\dot{\Phi}(t)=k_{B}\Bigg{[}\Gamma^{a+}_{\mathbb{AB}}\ln(\frac{{{\rm{k}}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}})+\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}}{{\rm{k}}_{\mathbb{BA}}^{b-}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}}{{\rm{k}}_{\mathbb{DB}}^{r-}})+\Gamma^{a-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{CD}}^{a+}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}}{{\rm{k}}_{\mathbb{CD}}^{b+}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}\\\
=&k_{B}\Bigg{[}\Gamma^{ab+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}}{{\rm{k}}_{\mathbb{BA}}^{a-}})+\Gamma^{r+}_{\mathbb{BD}}\ln(\frac{{\rm{k}}_{\mathbb{BD}}^{r+}}{{\rm{k}}_{\mathbb{DB}}^{r-}})+\Gamma^{ab-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{CD}}^{a+}})+\Gamma^{r-}_{\mathbb{CA}}\ln(\frac{{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}+k_{B}\Bigg{[}\Gamma^{b+}_{\mathbb{AB}}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}})\Bigg{]}\\\
=&k_{B}\Gamma_{\circlearrowright}\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}})+k_{B}(\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{b-}_{\mathbb{DC}})\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+k_{B}\Gamma^{b-}_{\mathbb{DC}}\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}})\\\
=&\kappa\Gamma_{\circlearrowright}\Bigg{[}\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{BD}}^{r+}{\rm{k}}_{\mathbb{DC}}^{a-}{\rm{k}}_{\mathbb{CA}}^{r-}}{{\rm{k}}_{\mathbb{BA}}^{a-}{\rm{k}}_{\mathbb{DB}}^{r-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AC}}^{r+}})\Bigg{]}+\\{\varepsilon_{\rm{L}}\Gamma^{b+}_{\mathbb{AB}}-(\varepsilon_{\rm{L}}+\kappa)\Gamma^{b-}_{\mathbb{DC}}\\}\Bigg{[}\left(\frac{k_{B}}{\kappa}\right)\ln(\frac{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}}{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}})\Bigg{]}\\\
+&(\Gamma^{b+}_{\mathbb{AB}}-\Gamma^{b-}_{\mathbb{DC}})\Bigg{[}k_{B}(1+\theta)\ln(\frac{{\rm{k}}_{\mathbb{AB}}^{b+}{\rm{k}}_{\mathbb{BA}}^{a-}}{{\rm{k}}_{\mathbb{BA}}^{b-}{\rm{k}}_{\mathbb{AB}}^{a+}})+k_{B}\theta\ln(\frac{{\rm{k}}_{\mathbb{DC}}^{b-}{\rm{k}}_{\mathbb{CD}}^{a+}}{{\rm{k}}_{\mathbb{CD}}^{b+}{\rm{k}}_{\mathbb{DC}}^{a-}})\Bigg{]}.\end{split}$
(94)