# Universal behaviour of Coulomb coupled Fermionic thermal diode

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

###### Abstract

We propose a minimal model of a Coulomb-coupled fermionic quantum dot thermal
diode that can act as an efficient thermal switch and exhibit complete
rectification behavior, even in the presence of a small temperature gradient.
Using two well-defined dimensionless system parameters, universal
characteristics of the optimal heat current conditions are identified. It is
shown to be independent of any system parameter and is obtained only at the
mean transitions point “$-0.5$”, associated with the equilibrium distribution
of the two fermionic reservoirs, tacitly referred to as “universal magic
mean”.

##  I Introduction

Heat management devices have attracted recent interest in nanoscale systems to
prevent overheating due to heat flow in desired areas of electronic circuits
Giazotto  _et al._ (2006); Pekola (2015); Benenti  _et al._ (2017); Roberts
and Walker (2011); Li  _et al._ (2012); Naseem  _et al._ (2020); Zeng and Wang
(2008); Ruokola  _et al._ (2009); Kuo and Chang (2010); Liu  _et al._ (2014);
Sánchez  _et al._ (2017); Landi  _et al._ (2021). Following the theoretical
proposal of a quantum dot thermal diode by Roukola and Ojanen Ruokola and
Ojanen (2011), a number of studies are carried out in achieving rectification
effect using simple quantum systems Terraneo  _et al._ (2002); Li  _et al._
(2004); Segal and Nitzan (2005); Ojanen (2009); Wu and Segal (2009); Werlang
_et al._ (2014); Mascarenhas  _et al._ (2014); Marcos-Vicioso  _et al._
(2018); Bhandari  _et al._ (2021); Iorio  _et al._ (2021); Balachandran  _et
al._ (2019); Kargı  _et al._ (2019); Ordonez-Miranda  _et al._ (2017); Iorio
_et al._ (2021); Aligia  _et al._ (2020); Upadhyay  _et al._ (2021); Díaz and
Sánchez (2021); Tesser  _et al._ (2022). Yet, most of the works that have been
done till now rely on either temperature gradient of bosonic reservoirs or
different coupling strengths between the system and the bath to break the
inversion symmetry of the overall system Tesser  _et al._ (2022). For example,
Werlang et. al. Werlang  _et al._ (2014) explored heat transport under the
influence of strong coupling between two spins interacting with their
respective bosonic bath. Miranda et. al. Ordonez-Miranda  _et al._ (2017)
identified similar diode characteristic in presence of different excitation
frequencies between the coupled spins. Based on a two and four terminal
quantum dot setups, Tesser et.al. Tesser  _et al._ (2022) recently explored
the role of level degeneracy and temperature bias. However, much less
attention is paid in achieving rectification effect by means of statistical
properties of the reservoir. Here, we provide a general framework to capture
the invariant aspect of the fermionic diodes in terms of two dimensionless
physical parameters and statistical distribution of the reservoir — the
uniqueness of which is found to be independent of the system energy levels,
interaction strength, bath spectrum, its temperature and chemical potentials.

To showcase our findings, we consider a Coulomb coupled quantum dot system and
study the interaction with two fermionic reservoirs with different temperature
and chemical potential. We find while temperature gradient governs the overall
heat flow direction, magnitude of heat current is primarily controlled by
chemical potential gradient. In contrast to bosonic counterpart, fermionic
rectifier allows us to control the heat current and switching effects in a
much more efficient way, even in presence of tiny temperature difference. Most
remarkably, we identify universal nature of complete, partial and no
rectification conditions that are valid for all Coulomb coupled fermionic
diodes.

The present work is organized as follows: We introduce model and dynamics in
Sec.II, and the steady state heat current in Sec. III. Microscopic picture
behind the thermodynamically consistent heat flow direction is summarized in
Sec. IV. While universal characteristics based on two dimensionless parameters
are presented in Sec. V, the role of efficient thermal switch and ideal
rectification effects are discussed in Sec. VI. Finally, we conclude in Sec.
VII.

##  II Model and Dynamics

Figure 1: Coulomb coupled QDs connected with fermionic reservoirs through
sequential tunnelling.

Our model consists of two quantum dots (QDs) which are strongly and
capacitively coupled to each other and interact through a long-range Coulomb
force so that they can exchange energy but no particles. We consider this
model following Ref Zhang  _et al._ (2017); Ruokola and Ojanen (2011); Zhang
_et al._ (2018); Tesser  _et al._ (2022); Aligia  _et al._ (2020) which have
been recently introduced to study the thermal diode and transistor effects in
a wide variety of Coulomb blockade quantum dot devices. We further assume that
each QD is tunnel-coupled to its fermionic reservoir Ghosh  _et al._ (2012).
While electron transport between QDs are forbidden due to Coulomb blockade
Zhang  _et al._ (2017), electron tunnelling between QDs and its respective
reservoir permits heat flow from one reservoir to another through the coupled
QD system. Within sequential tunnelling Zhang  _et al._ (2018) under the
Coulomb blockade regime, each QD can have only two-levels with occupation
number either zero or one. The two QDs, as well as the temperature and
chemical potential of the baths to which they are connected, are labelled by
indices $L$ and $R$ [Fig. 1]. The Hamiltonian of the coupled QDs is then given
by,

$\displaystyle
H_{D}=\sum_{\alpha=L,R}\varepsilon_{\alpha}d^{\dagger}_{\alpha}d_{\alpha}+\sum_{\alpha,\beta=L,R;\alpha\neq\beta}U_{\alpha\beta}d^{\dagger}_{\alpha}d_{\alpha}d^{\dagger}_{\beta}d_{\beta}.$
(1)

Here, $\varepsilon_{\alpha}$ is the lowest single-particle energy of QDs.
Without loss of any generality, we further assume
$\varepsilon_{L}<\varepsilon_{R}$. Here $U_{\alpha\beta}$ is the positive
Coulomb interaction energy between the electrons in different QDs, and
$d^{\dagger}_{\alpha}(d_{\alpha})$ denotes the creation (annihilation)
operator for the $\alpha$-th QD, whose eigenstates are $|0\rangle$ and
$|1\rangle$ with eigenvalues $0$ and $\varepsilon_{\alpha}$ respectively.
Since the interaction energy in Eq. (1) is diagonal in the eigenbasis of the
individual QDs, eigenstates of $H_{D}$ w.l.o.g can be written in terms of the
eigenstates of the two QDs in decreasing energy order, as
$|1\rangle=|00\rangle,|2\rangle=|10\rangle,|3\rangle=|01\rangle,|4\rangle=|11\rangle$.
For details we refer to Ref Zhang  _et al._ (2017); Ruokola and Ojanen (2011);
Zhang  _et al._ (2018). The Hamiltonian of the $\alpha$-th fermionic reservoir
is defined as $H_{R}^{\alpha}=\sum_{k}(\varepsilon_{k}-\mu_{\alpha})c_{\alpha
k}^{\dagger}c_{\alpha k}$ Tesser  _et al._ (2022), where $\varepsilon_{k}$ is
the energy of the non-interacting reservoir electrons with continuous wave
number $k$, $\mu_{\alpha}$ being the chemical potential and $c^{\dagger}(c)$
represents the creation (annihilation) operator of the electron reservoir. The
coupling between QDs and the respective reservoir is described by the
tunnelling Hamiltonian, $H_{T}^{\alpha}=\sum_{k}(t_{\alpha k}c_{\alpha
k}^{\dagger}d_{\alpha}+t_{\alpha k}^{*}d_{\alpha}^{\dagger}c_{\alpha k})$,
where $t_{\alpha k}$ is the tunnelling amplitude. Interaction within
sequential tunnelling approximation imposes restriction on simultaneous
tunnelling of more than one electrons at a time Zhang  _et al._ (2017); Tesser
_et al._ (2022); Ruokola and Ojanen (2011); Aligia  _et al._ (2020).
Consequently, there are in total four authorized transitions: the left
reservoir ($L$) induces transitions between $1\leftrightarrow 2$ and
$3\leftrightarrow 4$, while right bath ($R$) drives transitions between
$1\leftrightarrow 3$ and $2\leftrightarrow 4$. We define transition energies
$\omega_{ij}=\epsilon_{i}-\epsilon_{j}$, for $i>j$, where $\epsilon_{k}$ is
the eigenvalue of $H_{D}$ for the eigenstate $|k\rangle$. In the present case,
they read as $\omega_{21}=\varepsilon_{L}$, $\omega_{42}=\varepsilon_{R}+U$,
$\omega_{43}=\varepsilon_{L}+U$, and $\omega_{31}=\varepsilon_{R}$. The rates
at which the above transitions occur are computed using Lindblad master
equation Breuer and Petruccione (2002); Carmichael (2002).

We implement the strong-coupling formalism following Ref. Werlang  _et al._
(2014); Joulain  _et al._ (2016); Katz and Kosloff (2016); Jiang  _et al._
(2015); Goury and Sánchez (2019); Liu  _et al._ (2022); Gupt  _et al._ (2022)
to arrive at the master equation describing the time-evolution of the density
matrix $\rho$ of the coupled QDs system (See Appendix A)

$\frac{d\rho}{dt}=-\frac{i}{\hbar}[H_{D},\rho]+\mathcal{L}_{L}[\rho]+\mathcal{L}_{R}[\rho],$
(2)

under Born, Markov and secular approximation Werlang  _et al._ (2014); Marcos-
Vicioso  _et al._ (2018); Breuer and Petruccione (2002). It is important to
note here that the strong coupling formalism refers to the coupling between
the dots, while the the system-bath coupling is still assumed to be weak, so
that Born-Markov approximation can safely be implemented Werlang  _et al._
(2014); Joulain  _et al._ (2016); Katz and Kosloff (2016); Jiang  _et al._
(2015); Goury and Sánchez (2019); Liu  _et al._ (2022); Gupt  _et al._ (2022).
This implies that the Lindbladians are obtained on the basis of the
eigenstates of the full system Hamiltonian $H_{D}$. Thus dissipation mechanism
of each QD depends not only on the coupling to its own bath, but also on the
coupling between QDs, which is necessary for accurately describing the heat
flow and rectification effects over a wide range of system parameters as
considered below.

##  III Evaluation of Steady State Heat current

In the present model particles can not be exchanged between QDs as they
interact only through long range Coulomb force. As a result, there is no
particle flow in between reservoirs through the QD system Zhang  _et al._
(2017); Ruokola and Ojanen (2011); Zhang  _et al._ (2018); Tesser  _et al._
(2022); Aligia  _et al._ (2020). So, energy is exchanged only in the form of
heat [See Figs. 2 and 3]. The expression of the heat current can then be
obtained using the standard procedure Gelbwaser-Klimovsky  _et al._ (2015)
starting from the von-Neumann entropy of the system, defined as
$\mathcal{S}[\rho(t)]=-k_{B}{\rm Tr}[\rho(t)\ln\rho(t)]$. Upon taking the time
derivative of the von-Neumann entropy, one obtains

$\frac{d}{dt}\mathcal{S}[\rho(t)]=-\sum_{\alpha=L,R}k_{B}{\rm
Tr}\\{\mathcal{L}_{\alpha}[\rho(t)]\ln\rho(t)\\},$ (3)

where we have used the master equation (2) and ${\rm Tr}[\dot{\rho}]=0$. Sphon
inequality Spohn (1978) in the form of second law of thermodynamics Kosloff
(2013) for any Lindblad super-operator $\mathcal{L}$ can be written as ${\rm
Tr}\\{\mathcal{L}[\rho(t)](\ln\rho(t)-\ln\rho_{ss})\\}\leqslant 0$. Here
$\rho_{ss}$ is the steady state population of the system, satisfying
$\mathcal{L}[\rho_{ss}]=0$. Denoting the stationary state of
$\mathcal{L}_{\alpha}$, as $\rho_{ss}^{\alpha}$, above inequality can be
applied to both the terms of the sum in Eq. (3), yielding

$\sum_{\alpha=L,R}{\rm
Tr}\\{\mathcal{L}_{\alpha}[\rho(t)](\ln\rho(t)-\ln\rho_{ss}^{\alpha})\\}\leqslant
0.$ (4)

From Eqs. (3) and (4) we can write

$\begin{split}\frac{d}{dt}\mathcal{S}[\rho(t)]+\sum_{\alpha=L,R}k_{B}{\rm
Tr}[\mathcal{L}_{\alpha}[\rho(t)]\ln\rho_{ss}^{\alpha}]\geq 0\end{split}.$ (5)

Above equation can be compared with the dynamical version of the second law
given by Kosloff (2013); Gelbwaser-Klimovsky  _et al._ (2015)

$\frac{d}{dt}\mathcal{S}[\rho(t)]-\sum_{\alpha=L,R}\frac{J_{Q}^{\alpha}(t)}{T_{\alpha}}\geq
0.$ (6)

This allows us to identify the heat current (energy flow rate)
$J_{Q}^{\alpha}(t)$ associated with the $\alpha$-th reservoir as

$J_{Q}^{\alpha}(t)=-\frac{1}{\beta_{\alpha}}{\rm
Tr}\\{\mathcal{L}_{\alpha}[\rho(t)]\ln{\rho_{ss}^{\alpha}}\\}.$ (7)

For the steady state operation, the master equation (2) drives the system
towards a Gibbs-like stationary state, characterized by Gelbwaser-Klimovsky
_et al._ (2015); Gupt  _et al._ (2022)

$\rho_{ss}^{\alpha}=\mathcal{Z}^{-1}\exp[-\beta_{\alpha}H_{D}],$ (8)

where $\mathcal{Z}={\rm Tr}\\{\exp[-\beta_{\alpha}H_{D}]\\}$. For the detailed
derivation of the steady state heat current we refer to the review articles by
Kosloff and Gelbwaser et al. Kosloff (2013); Gelbwaser-Klimovsky  _et al._
(2015). Inserting Eq. (8) into Eq. (7), we obtain the general expression for
the heat current, which under stead state condition simplifies to

$J_{Q}^{\alpha}={\rm
Tr}\\{\mathcal{L}_{\alpha}[\rho_{ss}]H_{D}\\}=\sum_{\omega_{ij}}\omega_{ij}\Gamma^{\alpha}_{ij}.$
(9)

Here the net decaying rate $\Gamma_{ij}^{\alpha}$ from $|i\rangle$ to
$|j\rangle$ ($i>j$) is denoted as

$\begin{split}\Gamma_{ij}^{\alpha}=&-\Gamma_{ji}^{\alpha}\\\
\equiv&\gamma_{\alpha}[1-f_{\alpha}(\omega_{ij})]\rho_{ii}-\gamma_{\alpha}f_{\alpha}(\omega_{ij})\rho_{jj}=[\Gamma_{ij}^{\alpha}]_{\downarrow}-[\Gamma_{ij}^{\alpha}]_{\uparrow}.\end{split}$
(10)

The first term represents the emission and the second term corresponds to
absorption. $\gamma_{\alpha}$ being the bare tunnelling rate between dots and
respective reservoirs [Fig. 1]. $f_{\alpha}(\omega_{ij})$ is the Fermi
Distribution function (FDF)

$\begin{split}f_{\alpha}(\omega_{ij})\equiv
f_{\alpha}(\omega_{ij},T_{\alpha},\mu_{\alpha})=\bigg{[}{1+\exp\bigg({\frac{\omega_{ij}-\mu_{\alpha}}{k_{B}T_{\alpha}}})}\bigg{]}^{-1}\end{split},$
(11)

corresponding to the transition energy $\omega_{ij}=\epsilon_{i}-\epsilon_{j}$
between eigenstates $\rvert i\rangle$ and $\rvert j\rangle$ controlled by the
$\alpha$-th reservoir. Tunnelling of electrons into or out of QDs are
primarily governed by FDF. In our model there are four allowed transitions and
both lead guides two transitions each: $L$ lead drives transitions between
$\rvert 1\rangle\leftrightarrow\rvert 2\rangle$ and $\rvert
4\rangle\leftrightarrow\rvert 3\rangle$, while $R$ lead controls $\rvert
1\rangle\leftrightarrow\rvert 3\rangle$ and $\rvert
4\rangle\leftrightarrow\rvert 2\rangle$ transitions. For the sake of
convenience of our analysis, corresponding FDFs are expressed in terms of two
dimensionless parameters, effective tunnelling barrier
($\chi_{\alpha}={(\varepsilon_{\alpha}-\mu_{\alpha})}/{U}$) and dimensionless
thermal energy ($\xi_{\alpha}={k_{B}T_{\alpha}}/{U}$) as follows

$\begin{split}f_{L}(\omega_{21})=\left[{1+\exp\bigg({\frac{\varepsilon_{L}-\mu_{L}}{k_{B}T_{L}}}\bigg{missing})}\right]^{-1}=\left[{1+\exp\bigg({\frac{\chi_{L}}{\xi_{L}}}\bigg{missing})}\right]^{-1}=&f_{L}^{1},\\\
f_{L}(\omega_{43})=\left[{1+\exp\bigg({\frac{\varepsilon_{L}+U-\mu_{L}}{k_{B}T_{L}}}\bigg{missing})}\right]^{-1}=\left[{1+\exp\bigg({\frac{\chi_{L}+1}{\xi_{L}}}\bigg{missing})}\right]^{-1}=&f_{L}^{2},\\\
f_{R}(\omega_{31})=\left[{1+\exp\bigg({\frac{\varepsilon_{R}-\mu_{R}}{k_{B}T_{R}}}\bigg{missing})}\right]^{-1}=\left[{1+\exp\bigg({\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})}\right]^{-1}=&f_{R}^{1},\\\
f_{R}(\omega_{42})=\left[{1+\exp\bigg({\frac{\varepsilon_{R}+U-\mu_{R}}{k_{B}T_{R}}}\bigg{missing})}\right]^{-1}=\left[{1+\exp\bigg({\frac{\chi_{R}+1}{\xi_{R}}}\bigg{missing})}\right]^{-1}=&f_{R}^{2}.\end{split}$
(12)

Now, our task is to first evaluate the full expression of
$\Gamma_{ij}^{\alpha}$ at the steady state and then, find out the expression
of the heat current. Under the steady state condition master Eq. (2) is
characterized by $\dot{\rho}_{ss}=0$ which reduces to

$\begin{split}&\dot{\rho}_{11}=0=\Gamma_{31}^{R}-\Gamma_{12}^{L};\quad\quad\dot{\rho}_{22}=0=\Gamma_{12}^{L}-\Gamma_{24}^{R},\\\
&\dot{\rho}_{33}=0=\Gamma_{43}^{L}-\Gamma_{31}^{R};\quad\quad\dot{\rho}_{44}=0=\Gamma_{24}^{R}-\Gamma_{43}^{L}.\
\ \end{split}$ (13)

Thus, at the steady state, all the net transitions rates become equal to
$\Gamma$, i.e.,
$\Gamma_{31}^{R}=\Gamma_{12}^{L}=\Gamma_{24}^{R}=\Gamma_{43}^{L}\equiv\Gamma$.
The four set of equations in Eq. (LABEL:SS) are not independent since
$\sum_{i}\rho_{ii}=1$, which uniquely solves all the state occupation
probabilities as well as the heat current in terms of a single quantity
$\Gamma$:
$\Gamma_{31}^{R}=\Gamma_{12}^{L}=\Gamma_{24}^{R}=\Gamma_{43}^{L}\equiv\Gamma$.
From Eq. (9), we can then evaluate the general expression of heat current at
the steady state as follows

$\begin{split}J_{Q}^{R}=\varepsilon_{R}\Gamma_{13}^{R}-(\varepsilon_{R}+U)\Gamma_{42}^{R}=&-\varepsilon_{R}\Gamma+(\varepsilon_{R}+U)\Gamma=U\Gamma,\\\
J_{Q}^{L}=\varepsilon_{L}\Gamma_{12}^{L}-(\varepsilon_{L}+U)\Gamma_{43}^{L}=&\
\varepsilon_{L}\Gamma-(\varepsilon_{L}+U)\Gamma=-U\Gamma.\end{split}$ (14)

Thus we finally arrive at the explicit analytical expression of the steady
state heat current as

$\begin{split}J_{Q}^{L}=-J_{Q}^{R}=-U\Gamma.\end{split}$ (15)

To find out $\Gamma$, we rewrite Eq. (LABEL:SS) in terms of $f_{L(R)}^{1}$ and
$f_{L(R)}^{2}$

$\begin{split}\dot{\rho}_{11}=\gamma_{R}[1-f_{R}^{1}]\rho_{33}-[\gamma_{R}f_{R}^{1}+\gamma_{L}f_{L}^{1}]\rho_{11}+\gamma_{L}[1-f_{L}^{1}]\rho_{22}=0,&\\\
\dot{\rho}_{22}=\gamma_{L}f_{L}^{1}\rho_{11}-[\gamma_{L}(1-f_{L}^{1})+\gamma_{R}f_{R}^{2}]\rho_{22}+\gamma_{R}[1-f_{R}^{2}]\rho_{44}=0,&\\\
\dot{\rho}_{33}=\gamma_{R}f_{R}^{1}\rho_{11}-[\gamma_{L}f_{L}^{2}+\gamma_{R}(1-f_{R}^{1})]\rho_{33}+\gamma_{L}[1-f_{L}^{2}]\rho_{44}=0,&\\\
\dot{\rho}_{44}=\gamma_{L}f_{L}^{2}\rho_{33}-[\gamma_{R}(1-f_{R}^{2})+\gamma_{L}(1-f_{L}^{2})]\rho_{44}+\gamma_{R}f_{R}^{2}\rho_{22}=0.&\end{split}$
(16)

defined through Eq. (12) and find out the steady state populations subject to
the condition

$\rho_{11}+\rho_{22}+\rho_{33}+\rho_{44}=1.$ (17)

Using Eq. (16) and Eq. (17) we can construct

$\mathcal{M}\begin{bmatrix}\rho_{11}\\\ \rho_{22}\\\ \rho_{33}\\\
\rho_{44}\end{bmatrix}=\begin{bmatrix}0\\\ 0\\\ 0\\\ 1\end{bmatrix},$ (18)

where,

$\mathcal{M}=\begin{bmatrix}-[\gamma_{R}f_{R}^{1}+\gamma_{L}f_{L}^{1}]&\gamma_{L}[1-f_{L}^{1}]&\gamma_{R}[1-f_{R}^{1}]&0\\\
\gamma_{L}f_{L}^{1}&[\gamma_{L}(1-f_{L}^{1})+\gamma_{R}f_{R}^{2}]&0&\gamma_{R}[1-f_{R}^{2}]\\\
\gamma_{R}f_{R}^{1}&0&[\gamma_{L}f_{L}^{2}+\gamma_{R}(1-f_{R}^{1})]&\gamma_{L}[1-f_{L}^{2}]\\\
1&1&1&1\end{bmatrix}.$ (19)

Solving above equation, expressions for the steady state populations
$\rho_{ii}$ are found to be

$\begin{split}\rho_{11}=-\frac{1}{|\mathcal{M}|}[&\gamma_{L}{\gamma_{R}}^{2}\\{(1-f_{L}^{1})(1-f_{R}^{1})(1-f_{R}^{2})+(1-f_{L}^{2})(1-f_{R}^{1})f_{R}^{2}\\}\\\
+&\gamma_{R}{\gamma_{L}}^{2}\\{(1-f_{L}^{1})(1-f_{R}^{2})f_{L}^{2}+(1-f_{L}^{2})(1-f_{R}^{1})(1-f_{L}^{1})\\}],\\\
\rho_{22}=-\frac{1}{|\mathcal{M}|}[&\gamma_{L}{\gamma_{R}}^{2}\\{f_{L}^{1}(1-f_{R}^{1})(1-f_{R}^{2})+f_{L}^{2}(1-f_{R}^{2})f_{R}^{1}\\}\\\
+&\gamma_{R}{\gamma_{L}}^{2}\\{f_{L}^{1}(1-f_{R}^{2})f_{L}^{2}+(1-f_{L}^{2})(1-f_{R}^{1})f_{L}^{1}\\}],\\\
\rho_{33}=-\frac{1}{|\mathcal{M}|}[&\gamma_{L}{\gamma_{R}}^{2}\\{(1-f_{L}^{2})f_{R}^{1}f_{R}^{2}+(1-f_{L}^{1})(1-f_{R}^{2})f_{R}^{1}\\}\\\
+&\gamma_{R}{\gamma_{L}}^{2}\\{f_{R}^{2}(1-f_{L}^{2})f_{L}^{1}+(1-f_{L}^{2})f_{R}^{1}(1-f_{L}^{1})\\}],\\\
\rho_{44}=-\frac{1}{|\mathcal{M}|}[&\gamma_{L}{\gamma_{R}}^{2}\\{f_{L}^{1}f_{R}^{2}(1-f_{R}^{1})+f_{L}^{2}f_{R}^{1}f_{R}^{2}\\}\\\
+&\gamma_{R}{\gamma_{L}}^{2}\\{f_{L}^{1}f_{R}^{2}f_{L}^{2}+(1-f_{L}^{1})f_{L}^{2}f_{R}^{1}\\}],\end{split}$

where $|\mathcal{M}|$ stands for the determinant of the matrix $\mathcal{M}$.
Using the above expressions we can evaluate the final expression of $\Gamma$
as follows

$\Gamma=\frac{\gamma_{L}+\gamma_{R}}{\gamma_{L}\gamma_{R}}\left[\frac{f_{R}^{1}f_{R}^{2}f_{L}^{2}-f_{R}^{1}f_{R}^{2}f_{L}^{1}+f_{L}^{1}f_{L}^{2}f_{R}^{1}-f_{L}^{1}f_{L}^{2}f_{R}^{2}+f_{L}^{1}f_{R}^{2}-f_{L}^{2}f_{R}^{1}}{f_{L}^{1}f_{R}^{1}+f_{L}^{2}f_{R}^{2}-f_{L}^{1}f_{R}^{2}-f_{L}^{2}f_{R}^{1}-1}\right].$
(20)

Using Eq. (20), final expression of the steady state heat current in terms of
dimensionless parameters $\chi_{\alpha}$ and $\xi_{\alpha}$ can be written as

${J}^{L}_{Q}=-{J}^{R}_{Q}=-\frac{U(\gamma_{L}+\gamma_{R})}{\gamma_{L}\gamma_{R}}\left[\frac{\exp\bigg({\frac{1}{\xi_{L}}}\bigg{missing})-\exp\bigg({\frac{1}{\xi_{R}}}\bigg{missing})}{X}\right],$
(21)

where,

$\begin{split}X=&\bigg{[}\exp\bigg({\frac{1}{\xi_{L}}\bigg{}}\bigg{)})\left[2+\exp\bigg({-\frac{\chi_{R}}{\xi_{R}}\bigg{}}\bigg{)})+\exp\bigg({\frac{\chi_{L}}{\xi_{L}}\bigg{}}\bigg{)})+\exp\bigg({\frac{\chi_{L}}{\xi_{L}}-\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})\right]\\\
+&\exp\bigg({\frac{1}{\xi_{R}}}\bigg{missing})\left[2+\exp\bigg({-\frac{\chi_{L}}{\xi_{L}}}\bigg{missing})+\exp\bigg({\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})+\exp\bigg({\frac{\chi_{R}}{\xi_{R}}-\frac{\chi_{L}}{\xi_{L}}}\bigg{missing})\right]\\\
+&\exp\bigg({\frac{1}{\xi_{L}}+\frac{1}{\xi_{R}}}\bigg{missing})\left[\exp\bigg({\frac{\chi_{R}}{\xi_{R}}}+\exp{\frac{\chi_{L}}{\xi_{L}}}\bigg{missing})+\exp\bigg({\frac{\chi_{L}}{\xi_{L}}+\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})\right]\\\
+&\left[\exp\bigg({-\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})+\exp\bigg({-\frac{\chi_{L}}{\xi_{L}}}\bigg{missing})+\exp\bigg({-\frac{\chi_{L}}{\xi_{L}}-\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})\right]\bigg{]}.\end{split}$
(22)

Figure 2: Thermodynamically consistent and inconsistent heat flow directions:
Probable paths for transition cycle under (a) $T_{L}<T_{R}$ ; (b)
$T_{R}<T_{L}$. In both conditions, Path-I is initiated by the hot bath and
Path-II is initiated by the cold bath. Though Path-I seems to be natural heat
flow direction, transition cycle actually encompasses Path-II in both limits,
following the laws of thermodynamics.

This is the exact analytical expression of heat current derived under Born-
Markov master equation. Since $U$, $\gamma_{L}$, $\gamma_{R}$ and $X$ are all
positive, it is immediately clear from Eqs. (21)-(22), that if
$\xi_{R}>\xi_{L}$ ($T_{R}>T_{L}$), heat will flow from right to left
($J_{R}>0$) in accordance with Eq. (15).

##  IV Microscopic description of heat flow

Classically heat flows according to laws of thermodynamics from high to low
temperature. Quantum mechanically, flow of heat current must be governed by a
microscopic description, without violating the ultimate laws of thermodynamics
Kosloff (2013); Deffner and Campbell (2019); Ghosh  _et al._ (2019); Gupt  _et
al._ (2021); Sinha  _et al._ (2013). Since the QDs are strongly coupled with
each other, they form a four level system as depicted in Fig. 2, where inter-
dot transition is restricted due to long range Coulomb force. So, the energy
transfer is only caused by the heat baths via the coupled states of the
overall system Zhang  _et al._ (2017); Ruokola and Ojanen (2011); Zhang  _et
al._ (2018); Tesser  _et al._ (2022); Aligia  _et al._ (2020).For
$T_{R}>T_{L}$, natural heat flow direction will be from $R$ to $L$ with
$J_{Q}^{R}>0$. In this case, the first excitation from the ground state must
be guided by the cold bath, instead of hot bath, which may appear paradoxical
at first sight but makes the net transition rate $\Gamma>0$ as required by Eq.
(15). If the cycle follows the opposite path $1\rightarrow 3\rightarrow
4\rightarrow 2\rightarrow 1$ [Fig. 2(a)-I], then $J_{Q}^{L}>0$, as $U\Gamma$
amount of energy should be delivered by the left bath to the right bath, which
is in complete disagreement with the laws of thermodynamics. So, in order to
be consistent with the laws of the thermodynamics, transition cycle must run
in $1\rightarrow 2\rightarrow 4\rightarrow 3\rightarrow 1$, initiated by the
cold bath. This is seemingly paradoxical in the sense that classically we
expect during heat flow, energy is supplied by the hot bath and dumped into
the cold bath. For the specific example, it is however more favourable for the
cold bath to make the $1\rightarrow 2$ transition in Fig. 2(a)-II which costs
$\varepsilon_{L}$ amount of energy than $3\rightarrow 4$ transition in
Fig.2(a)-I which requires $(\varepsilon_{L}+U)$ amount of energy. It is
interesting to note that for $T_{L}>T_{R}$, (i.e., $J_{Q}^{L}>0$, heat flows
from left to right following Fig. 2(b)-II), first transition is still mediated
by the cold bath between $1\rightarrow 3$, as it requires less amount of
energy ($\varepsilon_{R}$) than $2\rightarrow 4$ transition
($\varepsilon_{R}+U$) in Fig. 2(b)-I.

Figure 3: Thermodynamically consistent transition cycle for $T_{R}>T_{L}$ in
case of $\varepsilon_{L}=\varepsilon_{R}=\varepsilon$; first excitation from
$|1\rangle$ must be guided by cold bath (Path-II: solid lines), instead of hot
bath (Path-I: dotted lines).

It is important to emphasize that although we have taken
$\varepsilon_{L}<\varepsilon_{R}$ to draw the schematic energy level diagram
of Fig. 2, magnitude of the heat current in particular depends only on the two
dimensionless parameters $\\{\chi_{\alpha},\xi_{\alpha}\\}$ via Eq. (21). In
the next section, we will explore the universal characteristics of the heat
current solely based on these two dimensionless system parameters. For
instance, Fig. 3 shows that basic principle behind the thermodynamically
consistent transition cycle remains same even in case of
$\varepsilon_{L}=\varepsilon_{R}$.

##  V Universal characteristic due to magic mean

While temperature gradient dictates the overall heat flow direction, chemical
potential plays a very important role in determining the heat current
magnitudes. To illustrate, we note that FDFs $f_{L(R)}^{1}$ and $f_{L(R)}^{2}$
in Eq. (12), expressed in terms of dimensionless thermal energy
$\xi_{L(R)}=k_{B}T_{L(R)}/U$ and effective tunnelling barrier
$\chi_{L(R)}=(\varepsilon_{L(R)}-\mu_{\alpha})/U$, are constrained by $0\leq
f_{L(R)}^{2}<f_{L(R)}^{1}\leq 1$ and become $0.5$ at $\chi_{L(R)}=0$ and $-1$,
which correspond to $\mu_{L(R)}=\omega_{21(31)}=\varepsilon_{L(R)}$ and
$\mu_{L(R)}=\omega_{43(42)}=\varepsilon_{L(R)}+U$, respectively. According to
Eq. (10), if $0\leq\\{f_{L(R)}^{1},f_{L(R)}^{2}\\}\leq 0.5$, it favours de-
excitation and equivalently excitation is favoured when
$0.5\leq\\{f_{L(R)}^{1},f_{L(R)}^{2}\\}\leq 1$, for the corresponding
transition. This is unique to fermionic reservoirs and allows us to implement
a systematic analytical scheme solely based on FDFs which is in sharp contrast
to its bosonic counterpart. The span of $f_{L(R)}^{1}$ and $f_{L(R)}^{2}$ as a
function $\chi_{L(R)}$, can thus be classified into five domains
$\mathcal{D}_{[1-5]}$ [Fig. 4: Main] and the total FDF, defined as
$f_{L(R)}=f_{L(R)}^{1}+f_{L(R)}^{2}$, is found to be spreading between $0\leq
f_{L(R)}\leq 2$ [Fig. 4: Inset].

Figure 4: $\mathcal{D}_{3}$ always spreads between transition points of
$f_{L(R)}^{1}$ and $f_{L(R)}^{2}$; span of $\mathcal{D}_{1[5]}$ and
$\mathcal{D}_{2[4]}$ are not fixed and strongly depends on the value of
$\xi_{L(R)}$. Color gradients signifies the magnitude of $|J_{Q}|$.

In domain $\mathcal{D}_{1[5]}$, both $f_{L(R)}^{1}$ and $f_{L(R)}^{2}\sim
1[0]$, thus Eq. (10) reduces to
$\Gamma_{21/43}^{L}\simeq\\{\Gamma_{21/43}^{L}\\}_{\uparrow[\downarrow]}$ and
$\Gamma_{31/42}^{R}\simeq\\{\Gamma_{31/42}^{R}\\}_{\uparrow[\downarrow]}$ i.e.
only absorption (emission) is allowed in all four transitions, hence
transition cycle can't be completed [Fig. 2] which results in vanishing
$\rvert J_{Q}\arrowvert\rightarrow 0$. Now, $\mathcal{D}_{2[4]}$ is
characterised by $0.5[0]\lesssim\\{f_{L(R)}^{1},f_{L(R)}^{2}\\}\lesssim
1[0.5]$, hence excitation as well as de-excitation is favoured for all four
transitions yielding non-zero heat current. On the contrary, domain
$\mathcal{D}_{3}$, parametrized by $-1\leq\chi_{L(R)}\leq 0$, is sharply
defined between the transition points $f_{L(R)}^{2}$ and $f_{L(R)}^{1}$. As
opposed to $\mathcal{D}_{3}$, both $f_{L(R)}^{1}$ and $f_{L(R)}^{2}$ are
closer to $1[0]$ in $\mathcal{D}_{2[4]}$, so that heat flux decreases in
$\mathcal{D}_{2[4]}$ relative to $\mathcal{D}_{3}$. In $\mathcal{D}_{3}$, both
$f_{L(R)}^{1}$ and $f_{L(R)}^{2}$ take values such that excitation from
$\lvert 1\rangle\rightarrow\lvert 2\rangle(\lvert 2\rangle\rightarrow\lvert
4\rangle)$ and de-excitation from $\lvert 4\rangle\rightarrow\lvert
3\rangle(\lvert 3\rangle\rightarrow\lvert 1\rangle)$ can occur simultaneously
at an optimal rate. Condition becomes ideal at the midpoint of
$\mathcal{D}_{3}$ (point $B$ in Fig. 4) which corresponds to
$\chi_{L(R)}=-0.5$. To find out analytically, the optimal value of
$\chi_{L(R)}$ for which $\arrowvert{J_{Q}}\rvert$ becomes maximum, first we
have to differentiate ${J_{Q}}\equiv J_{Q}^{R}$ (Cf. Eq. (21)) w.r.t.
$\chi_{L(R)}$ and set that equal to be zero. Now, by differentiating Eq. (21)
w.r.t $\chi_{R}$, for fixed $\\{\chi_{L},\xi_{L},\xi_{R}\\}$, we obtain,

$\begin{split}&\left(\frac{\partial
J_{Q}^{R}}{\partial\chi_{R}}\right)_{\chi_{L},\xi_{L},\xi_{R}}\\\
&=\frac{U(\gamma_{L}+\gamma_{R})}{\gamma_{L}\gamma_{R}}\left[\frac{\partial}{\partial\chi_{R}}\left\\{\frac{f_{R}^{1}f_{R}^{2}f_{L}^{2}-f_{R}^{1}f_{R}^{2}f_{L}^{1}+f_{L}^{1}f_{L}^{2}f_{R}^{1}-f_{L}^{1}f_{L}^{2}f_{R}^{2}+f_{L}^{1}f_{R}^{2}-f_{L}^{2}f_{R}^{1}}{f_{L}^{1}f_{R}^{1}+f_{L}^{2}f_{R}^{2}-f_{L}^{1}f_{R}^{2}-f_{L}^{2}f_{R}^{1}-1}\right\\}\right]\\\
&=0.\end{split}$ (23)

As, $U>0$ and $\gamma_{L},\gamma_{R}\neq 0$, Eq. (LABEL:A17) implies that the
only non-trivial solution for maximizing $|J_{Q}|$ is equivalent to maximizing
$\Gamma$, which is given by the criteria

$f_{R}^{1}+f_{R}^{2}=1.$ (24)

It is clear from Eq. (LABEL:A17) that heat current vanishes under two limiting
conditions: $(i)\ f_{R}^{1}=f_{R}^{2}=0;\;(ii)\ f_{R}^{1}=f_{R}^{2}=1$. Again
from Eq. (12), we can write

$1\geq f_{\alpha}^{1}\geq f_{\alpha}^{2}\geq 0,\quad\alpha=L,R.$ (25)

So, Eq. (25) signifies that, if $f_{\alpha}^{1}=0$ then $f_{\alpha}^{2}$ is
certainly $0$. Similarly, if $f_{\alpha}^{2}=1$ then $f_{\alpha}^{1}$ is
certainly $1$. In both cases, heat current vanishes. Thus, the only condition
for nonzero heat current reduces to $f_{R}^{1}\neq 0$ and $f_{R}^{2}\neq 1$.
Under these conditions, after solving Eq. (24), we get the value of
$\chi_{R}=-0.5$, for which $|J_{Q}|$ attains maximum:

$\begin{split}\ &f_{R}^{1}=1-f_{R}^{2}\\\ \text{or},\
&\frac{1}{f_{R}^{1}}=\frac{1}{1-f_{R}^{2}}\\\ \text{or},\
&1+\exp\bigg({\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})=1+\exp\bigg({-{\frac{\chi_{R}+1}{\xi_{R}}}}\bigg{missing})\\\
\text{or},\
&\exp\bigg({\frac{\chi_{R}}{\xi_{R}}}\bigg{missing})=\exp\bigg({-{\frac{\chi_{R}+1}{\xi_{R}}}}\bigg{missing})\\\
\text{or},\ &\frac{\chi_{R}}{\xi_{R}}=-{\frac{\chi_{R}+1}{\xi_{R}}}\\\
\text{or},\ &\chi_{R}=-\chi_{R}-1\\\ \text{or},\ &\chi_{R}=-0.5.\end{split}$
(26)

So, heat current will be maximum when $\chi_{R}=-0.5$ [Fig 5(a)]; similarly,
when $\chi_{L}$ varies, $\chi_{L}=-0.5$ is the criteria for having highest
heat current. Hence, we can conclude that $\rvert J_{Q}\arrowvert$ is maximum
when both $\chi_{L(R)}=-0.5$, irrespective of $\xi_{L(R)}$, which is supported
by numerical simulations [Fig. 5(b,c)] and also follows from analytically
derived condition $f_{L(R)}=1$ [Cf. Eq. (24)] for maximum heat current [Fig.
4:Inset]. Now, putting these conditions in Eq. (21), we can evaluate the exact
analytical expression of the maximum heat current as

$\begin{split}&|J_{Q}|_{max}\\\
&=\frac{U(\gamma_{L}+\gamma_{R})}{\gamma_{L}\gamma_{R}}\bigg{\rvert}\frac{\sinh{(\frac{1}{2\xi_{L}}-\frac{1}{2\xi_{R}})}}{2[1+\cosh(\frac{1}{2\xi_{R}})+\cosh(\frac{1}{2\xi_{L}})]+\cosh{(\frac{1}{2\xi_{L}}-\frac{1}{2\xi_{R}})}+\cosh{(\frac{1}{2\xi_{R}}-\frac{1}{2\xi_{L}})}}\bigg{\arrowvert}\\\
&\equiv U|\Gamma|_{max}.\end{split}$ (27)

Figure 5: (a) Variation of heat current $J_{Q}^{R}$ (red line) with $\chi_{R}$
for fixed $\chi_{L}$, $\xi_{R}$ and $\xi_{L}$. Heat current become maximum at
$\chi_{R}=-0.5$ which is supported by the plots of first and second
derivatives of heat current: $J_{Q}^{\prime R}$ (green line) is zero and
$J_{Q}^{\prime\prime R}$ (blue line) is negative at the magic mean point
$\chi_{R}=-0.5$. Absolute value of the heat current $\arrowvert J_{Q}\rvert$
is plotted as a function of (b) $\\{\chi_{L},\chi_{R}\\}$ for $\xi_{R}=0.5$
and $\xi_{L}=0.2$; (c) $\\{\chi_{L},\chi_{R}\\}$ for $\xi_{R}=1$ and
$\xi_{L}=0.5$

Several remarks are now in order:

  * •

Maximum heat current is obtained at $\chi_{L(R)}=-0.5$ and the magnitude only
depends on the temperature of two leads, their tunnelling rates, and the
Coulomb interaction between the dots. Since, $-0.5$ is the mean of $0$ and
$-1$ which are transition points of $f_{L(R)}^{1}$ and $f_{L(R)}^{2}$
respectively (point $B$ in Fig. 4:main) and this mean is also independent of
other controlling parameters, we term this number “$-0.5$” as the “universal
magic mean”. Significance of the point $B$ lies on the fact that at the magic
mean $\chi_{L(R)}=-0.5$, chemical potential of left (right) lead becomes
exactly equal to mean of the transition energies $\omega_{21(31)}$ and
$\omega_{43(42)}$ driven by that bath:

$\begin{split}&\chi_{L(R)}=-0.5\\\
\Rightarrow\quad&\varepsilon_{L(R)}-\mu_{L(R)}=-0.5U\\\
\Rightarrow\quad&\mu_{L(R)}=\frac{2\varepsilon_{L(R)}+U}{2}\\\
\Rightarrow\quad&\mu_{L}=\frac{\omega_{43}+\omega_{21}}{2};\quad\mu_{R}=\frac{\omega_{42}+\omega_{31}}{2}.\end{split}$
(28)

and therefore provides maximum control over $L(R)$ bath to guide both the
absorption and decay simultaneously at maximum $\Gamma$, resulting maximum
$\rvert J_{Q}\arrowvert$.

  * •

With the increase of $\xi_{L(R)}$, $|J_{Q}|$ also spreads out, keeping maxima
point fixed at $\chi_{L(R)}=-0.5$. This precisely indicates that maximum heat
current will invariably be obtained at $\chi_{L(R)}=-0.5$, irrespective of all
$\xi_{L,R}$ [Fig. 5(b,c)]. 3D plots along $\chi_{L}$ are squeezed for smaller
values of $\xi_{L(R)}$[Fig. 5(b)] (Assuming $\xi_{L}<\xi_{R}$ as in Fig. 2)
and if we increase $\xi_{L}(\xi_{R})$, it gets expanded along both positive
and negative $\chi_{L(R)}$, leaving out point of maxima intact at
$\chi_{L(R)}=-0.5$ [Fig. 5(c)]. This fundamental feature of the “magic mean”
makes it truly universal.

##  VI Efficient Heat current modulator

With the increase of $\xi_{L(R)}$, $\mathcal{D}_{1[5]}$ gets narrower while
$\mathcal{D}_{2[4]}$ spreads out, without affecting $\mathcal{D}_{3}$ [Fig.
4]. As a consequence, $\arrowvert J_{Q}\rvert$ spreads out with $\xi_{L(R)}$,
keeping point of maxima fixed at $\chi_{L(R)}=-0.5$. Now if $\xi_{L(R)}\ll 1$,
then the span of $\mathcal{D}_{2[4]}$ is very small compared to
$\mathcal{D}_{1[5]}$ while domain $\mathcal{D}_{3}$ always remains in between
$-1\leq\chi_{L(R)}\leq 0$. So, there are effectively two domains: (I) Domain
ON, where $\arrowvert J_{Q}\rvert\neq 0$ and (II) Domain OFF, where
$\arrowvert J_{Q}\rvert=0$ and switching effect gets more prominent if
$\xi_{L(R)}\ll 1$. Thus we can operate our model as an efficient thermal
switch to on-off heat current just by shifting the domain from ON to OFF
through the change of $\chi_{L(R)}$ or in turn controllable experimental
parameter $\varepsilon_{L(R)}$ [Fig. 6a]. This is the basic underlying
principle behind switching effect for all coulomb coupled fermionic thermal
diodes. It becomes more prominent for smaller cold bath temperature
$\xi=k_{B}T/U$ and can be achieved more easily by varying $U$, instead of
lowering $T$ to smaller value.

Figure 6: (a) Main: Thermal switching effect is obtained with the variation of
scaled heat current $|J_{Q}|/\gamma U$ w.r.t $\chi_{cold}$ and gets more
prominent with $\xi_{cold}\ll 1$ as Domain ON and OFF become more precise;
Inset: Expressing switching effect in terms of controllable parameter
$\varepsilon_{L(R)}$. (b) No rectification ($\mathcal{R}=0$) is obtained
irrespective of the value of $\Delta\xi$, if $\bar{\chi}=-0.5$ or
$\Delta{\chi}=0$ (red line). Partial (purple and green lines) and complete
(blue line) rectifications are obtained for nonzero $|\Delta\chi|$ and
$\bar{\chi}$ deviating from the magic mean $-0.5$. (c) Main: $\xi_{cold}$
(orange and green lines) has stronger dependence on rectification than
$\xi_{hot}$ (blue dashed and dashed-dot lines) for an equal change in
$|\Delta\xi|$, where solid blue line corresponds to the blue curve of Fig.6b;
Inset: Magnitude of heat current decreases by an amount $10^{-6}$ on reversing
the temperature gradient for the blue curve in Fig.6b.

Finally, the model operates as an efficient thermal diode with rectification
factor ($\mathcal{R}$) approaching $1$, even in presence of low temperature
gradient. This produces a clear advantage over all previously proposed models
with $\mathcal{R}$ factor defined as Bhandari  _et al._ (2021); Roberts and
Walker (2011); Segal and Nitzan (2005); Wu and Segal (2009); Joulain  _et al._
(2016),

$\mathcal{R}(\Delta\xi)=\frac{\arrowvert
J_{Q}^{R}(\Delta\xi)-J_{Q}^{L}(-\Delta\xi)\arrowvert}{\arrowvert
J_{Q}^{R}(\Delta\xi)+J_{Q}^{L}(-\Delta\xi)\lvert},$ (29)

where $\Delta\xi=\xi_{R}-\xi_{L}$ is identified as the temperature gradient,
for $\xi_{R}>\xi_{L}$. With $\Delta\xi>0$, i.e.
$\xi_{R(L)}\equiv\xi_{hot(cold)}$, heat current flows from $R$ to $L$
fulfilling $J_{Q}^{R}>0$ and if we exchange temperatures of the bath, then
temperature gradient becomes $-\Delta\xi$ or $\xi_{R(L)}\equiv\xi_{cold(hot)}$
and consequently heat flows from $L$ to $R$ satisfying $J_{Q}^{L}>0$. Now if
the heat current vanishes on reversing the temperature gradient, i.e.,
$\arrowvert J_{Q}\rvert$ is finite in one direction but null in the other,
then complete rectification is achieved with $\mathcal{R}\rightarrow 1$,
irrespective of the value of $\Delta\xi$, whereas $\mathcal{R}\rightarrow 0$
corresponds to no rectification i.e. no change in $\arrowvert J_{Q}\rvert$
upon inverting the temperature gradient. So, current asymmetry in two
directions is primary criteria behind positive rectification: (i) It is clear
from Fig. 4 that for a given $\Delta\xi$, $\arrowvert J_{Q}\rvert$ depends on
$|\Delta\chi|$. If $\chi_{hot}=\chi_{cold}$, or
$\Delta\chi=\chi_{hot}-\chi_{cold}=0$, heat currents will be symmetric in both
ways, therefore $\mathcal{R}=0$. (ii) From Fig. 6a, we find the variation of
$\arrowvert J_{Q}\rvert$ w.r.t $\chi_{\alpha}$, is completely symmetric about
the magic mean $-0.5$, i.e., $\rvert J_{Q}(\chi_{\alpha})\lvert=\rvert
J_{Q}(-1-\chi_{\alpha})\arrowvert$. As a result, following Eq. (29) and Fig.
6a, we finds if $\chi_{hot}=-1-\chi_{cold}$ or
$\bar{\chi}=\frac{1}{2}(\chi_{hot}+\chi_{cold})=-0.5$, then
$J_{Q}^{R}(\Delta\xi)=J_{Q}^{L}(-\Delta\xi)$, yielding $\mathcal{R}=0$. Thus,
$\mathcal{R}$ can be zero either (i) the difference between effective
tunnelling barriers $\Delta\chi=0$, or (ii) mean effective tunnelling barrier
is equal to magic mean ($\bar{\chi}=-0.5$); else rectification occurs.

Once $|\Delta\chi|\neq 0$ and $\bar{\chi}\neq-0.5$, effective tunnelling
barriers $\chi_{L(R)}$ make dissimilar effects on $|J_{Q}|$ upon reversing the
temperature gradient and consequently $\mathcal{R}$ becomes nonzero. For a
fixed value of $|\Delta\xi|$, $\mathcal{R}$ increases as $\Delta\chi$ shifts
from zero and approaches one as $\bar{\chi}$ significantly deviates from the
magic mean [Fig. 6b]. The reason behind this, for smaller $|\Delta\chi|$, the
effect of tunnelling barriers on $|J_{Q}|$ are comparable on reversing the
temperature gradient, yielding partial rectification ($0<\mathcal{R}<1$). But,
with the increase in $|\Delta\chi|$, the effect of $\chi_{L(R)}$ are no longer
compatible and it creates larger heat current asymmetry between the two
directions, leading to complete rectification behaviour with
$\mathcal{R}\rightarrow 1$. As heat current is symmetric about
$\bar{\chi}=-0.5$, we only consider positive variation of $\bar{\chi}$ in
Fig.6b. Moreover, variation of $\mathcal{R}$ with $|\Delta\chi|$ depends on
the value of $|\Delta\xi|$, larger $\Delta\xi$ implies better rectification.
As $\xi_{cold}$ is less than $\xi_{hot}$ and the transition cycle is always
initiated by the cold bath, heat current executes stronger dependence on
$\xi_{cold}$ than $\xi_{hot}$ for a given change in $|\Delta\xi|$. Thus the
variation of $\mathcal{R}$ with $|\Delta\chi|$ yields appreciable change when
a given $|\Delta\xi|$ is altered due to change in $\xi_{cold}$ than
$\xi_{hot}$ [Fig.6c: Main]. So, the complete rectification is more favourable
when $\xi_{cold}\ll 1$ [Fig. 6a]. We may have complete rectification even
without considering absolute temperature of the cold bath close to zero, as we
can vary $\xi_{cold}=k_{B}T/U$ by changing $U$ and keeping $T$ finite.
Finally, Fig. 6c[Inset] shows that heat currents become $\sim 10^{-6}$ times
smaller on reversing temperature gradient which corresponds to
$\mathcal{R}\approx 1$ curve in Fig.6b.

##  VII Conclusion

To conclude, the present model can be implemented as an efficient thermal
switch as well as thermal rectifier to modulate heat current close to ideal
rectification, even at arbitrarily low temperature differences between the
heat reservoirs. Independent of the details of the system, the heat flow and
heat rectification are characterized by a small set of universal parameters.
The position of maximum heat current at the magic mean “$-0.5$” in terms of
dimensionless physical parameters, is the major findings of the domain
analysis scheme presented here. The magic mean is robust and universal in the
sense that it is invariant w.r.t the variation of any other system or bath
parameters and truly reflects the impact of chemical potential to decide the
magnitudes of heat current. Present protocol is unique to fermionic systems
and can be applied to more complicated three or multi-terminal fermionic
devices. The straightforward generalisation of the present scheme to multi-
terminal set-up would be computing the magic mean associated with the multiple
fermionic reservoirs. The advantage of using quantum dot systems is that they
have discrete energy levels with strong on-site Coulomb interaction, and can
be simultaneously tunnel coupled to their respective reservoirs. Their
discrete energy levels provides energy-selective transport and can be tunned
via the application of the external gate voltages. In view of recent
experimental advances in Coulomb coupled quantum-dot systems Thierschmann  _et
al._ (2015a); Hartmann  _et al._ (2015); Thierschmann  _et al._ (2015b);
Pfeffer  _et al._ (2015), our findings will have important implications in
designing novel thermal devices and opening up potential applications in
controlling thermal current at nanoscales.

## Acknowledgements

This research was funded by the Initiation grant of IITK (Grant No.
IITK/CHM/2018513) and SRG, SERB (grant no. SRG/2019/000289) India. S.G. and
N.G. are thankful to CSIR for the fellowship. S.G. is also grateful to the
Ministry of Education, Government of India for the Prime Minister Research
Fellowship (PMRF).

##  Appendix A Derivation of the Lindblad Master equation

First we derive the master equation for our composite quantum dot system which
is coupled with fermionic reservoirs. Let us start with the total interaction
Hamiltonian given by

$H_{T}=\sum_{\alpha=L,R}H_{T}^{\alpha}=\sum_{\alpha=L,R}\sum_{k}(t_{\alpha
k}c_{\alpha k}^{\dagger}d_{\alpha}+t_{\alpha
k}^{*}d_{\alpha}^{\dagger}c_{\alpha k}).$ (30)

It should be noted that one can not write the fermionic interaction
Hamiltonian using a tensor product representation, since operators involved in
the tensor product commute by construction while fermionic operators
anticommute. To formally derive the master equation using tensor product
representation, one needs more formal approach in terms of Jordon-Wigner
transformation. For details we refer to Ref Schaller  _et al._ (2009); Zedler
_et al._ (2009). In the present case, starting from Eq. (30), we can derive
the von-Neumann equation for the total density operator of the composite
system $\rho_{T}(t)$, as

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

Integrating the above equation following Ref. Zedler  _et al._ (2009), one
obtains

$\begin{split}\frac{\partial}{\partial
t}\rho_{D}(t)=\frac{1}{(i\hbar)^{2}}\int_{0}^{t}dt^{\prime}{\rm
Tr}_{L,R}[H_{T}(t),[H_{T}(t^{\prime}),\rho_{T}(t^{\prime})]],\end{split}$ (32)

where we denote ${\rm Tr}_{L,R}\\{\rho_{T}(t)\\}=\rho_{D}(t)$ as the reduced
density operator for the system and also assumed that ${\rm
Tr}_{L,R}[H_{T}(t),\rho_{T}(0)]=0$. Here ${\rm Tr}_{L,R}$ refers to the trace
over each bath degrees of freedom. The reduced dynamics of the system in the
weak sequential tunnelling limit can then be written as Marcos-Vicioso  _et
al._ (2018); Wijesekara  _et al._ (2021); Carmichael (2002); Breuer and
Petruccione (2002)

$\begin{split}\dot{\rho}_{D}(t)=\frac{1}{(i\hbar)^{2}}\int_{0}^{\infty}dt^{\prime}{\rm
Tr}_{L,R}[H_{T}(t),[H_{T}(t-t^{\prime}),\rho_{D}(t)\otimes\rho_{L}\otimes\rho_{R}]],\end{split}$
(33)

where we substitute $\rho_{T}(t)=\rho_{D}(t)\otimes\rho_{L}\otimes\rho_{R}$.
Since the bath operator obeys ${\rm
Tr}_{\alpha}\\{c_{\alpha}(t)\rho_{\alpha}\\}=0={\rm
Tr}_{\alpha}\\{c^{{\dagger}}_{\alpha}(t)\rho_{\alpha}\\}$ for $\alpha=L,R$, we
obtain Breuer and Petruccione (2002); Aligia  _et al._ (2020)

${\rm
Tr}_{L,R}\\{[H_{T}^{\alpha}(t),[H_{T}^{\beta}(t-t^{\prime}),\rho_{D}(t)\otimes\rho_{L}\otimes\rho_{R}]]\\}=0\quad\quad\alpha\neq\beta;\alpha,\beta=L,R.$
(34)

As a result Eq. (33) simplifies to

$\dot{\rho}_{D}(t)=\frac{1}{(i\hbar)^{2}}\sum_{\alpha=L,R}\left\\{\int_{0}^{\infty}dt^{\prime}\
{\rm
Tr}_{L,R}[H_{T}^{\alpha}(t),[H_{T}^{\alpha}(t-t^{\prime}),\rho_{D}(t)\otimes\rho_{L}\otimes\rho_{R}]]\right\\}.$
(35)

Following the standard procedure of Zedler  _et al._ (2009); Ghosh  _et al._
(2012), one can then derive the master equation

$\dot{\rho}_{D}(t)=\mathcal{L}_{L}[\rho_{D}(t)]+\mathcal{L}_{R}[\rho_{D}(t)],$
(36)

where the Lindblad operators $\mathcal{L}_{\alpha}[\rho]$ are given by

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

The operator $d_{\alpha}(\omega_{\alpha})$ assumes the form of
$|i\rangle\langle j|$ ($i\neq j;\;i,j=1,2,3,4$) and causes the transition
driven by left (right) reservoir with positive energy $\omega_{\alpha}$ such
that $\omega_{L}=\omega_{21},\;\omega_{43}$ and
$\omega_{R}=\omega_{31},\;\omega_{42}$. In above Eq.(A) the temperature
dependent bath auto-correlations functions are given by Zedler  _et al._
(2009); Schaller  _et al._ (2009); Ghosh  _et al._ (2012),

$\displaystyle\mathcal{G}_{\alpha}(\omega_{\alpha})=$
$\displaystyle\frac{\gamma_{\alpha}}{2}f_{\alpha}(\omega_{\alpha});\quad\mathcal{G}_{\alpha}(-\omega_{\alpha})=$
$\displaystyle\frac{\gamma_{\alpha}}{2}(1-f_{\alpha}(\omega_{\alpha})).$ (38)

So, using the above relations in Eq. (A), we obtain the final expression of
the Lindbladian operator as

$\displaystyle\mathcal{L}_{\alpha}[\rho]$ $\displaystyle=$
$\displaystyle\sum_{\omega_{\alpha}>0}\frac{\gamma_{\alpha}}{2}f_{\alpha}(\omega_{\alpha})\left[d_{\alpha}^{\dagger}(\omega_{\alpha})\rho
d_{\alpha}(\omega_{\alpha})-\frac{1}{2}{\\{\rho},d_{\alpha}(\omega_{\alpha})d_{\alpha}^{\dagger}(\omega_{\alpha})\\}\right]$
$\displaystyle+$
$\displaystyle\frac{\gamma_{\alpha}}{2}(1-f_{\alpha}(\omega_{\alpha}))\left[d_{\alpha}(\omega_{\alpha})\
\rho\
d_{\alpha}^{\dagger}(\omega_{\alpha})-\frac{1}{2}{\\{\rho},d_{\alpha}^{\dagger}(\omega_{\alpha})d_{\alpha}(\omega_{\alpha})\\}\right].$

## References

  * Giazotto  _et al._ (2006) F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, Rev. Mod. Phys. 78, 217 (2006).
  * Pekola (2015) J. P. Pekola, Nature Physics 11, 118 (2015).
  * Benenti  _et al._ (2017) G. Benenti, G. Casati, K. Saito, and R. Whitney, Physics Reports 694, 1 (2017).
  * Roberts and Walker (2011) N. Roberts and D. Walker, International Journal of Thermal Sciences 50, 648 (2011).
  * Li  _et al._ (2012) N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * Naseem  _et al._ (2020) M. T. Naseem, A. Misra, O. E. Müstecaplioğlu, and G. Kurizki, Phys. Rev. Research 2, 033285 (2020).
  * Zeng and Wang (2008) N. Zeng and J.-S. Wang, Phys. Rev. B 78, 024305 (2008).
  * Ruokola  _et al._ (2009) T. Ruokola, T. Ojanen, and A.-P. Jauho, Phys. Rev. B 79, 144306 (2009).
  * Kuo and Chang (2010) D. M.-T. Kuo and Y.-c. Chang, Phys. Rev. B 81, 205321 (2010).
  * Liu  _et al._ (2014) Y.-Y. Liu, W.-X. Zhou, L.-M. Tang, and K.-Q. Chen, Applied Physics Letters 105, 203111 (2014).
  * Sánchez  _et al._ (2017) R. Sánchez, H. Thierschmann, and L. W. Molenkamp, New Journal of Physics 19, 113040 (2017).
  * Landi  _et al._ (2021) G. T. Landi, D. Poletti, and G. Schaller,  (2021), 10.48550/arxiv.2104.14350.
  * Ruokola and Ojanen (2011) T. Ruokola and T. Ojanen, Phys. Rev. B 83, 241404 (2011).
  * Terraneo  _et al._ (2002) M. Terraneo, M. Peyrard, and G. Casati, Phys. Rev. Lett. 88, 094302 (2002).
  * Li  _et al._ (2004) B. Li, L. Wang, and G. Casati, Phys. Rev. Lett. 93, 184301 (2004).
  * Segal and Nitzan (2005) D. Segal and A. Nitzan, Phys. Rev. Lett. 94, 034301 (2005).
  * Ojanen (2009) T. Ojanen, Phys. Rev. B 80, 180301 (2009).
  * Wu and Segal (2009) L.-A. Wu and D. Segal, Phys. Rev. Lett. 102, 095503 (2009).
  * Werlang  _et al._ (2014) T. Werlang, M. A. Marchiori, M. F. Cornelio, and D. Valente, Phys. Rev. E 89, 062109 (2014).
  * Mascarenhas  _et al._ (2014) E. Mascarenhas, D. Gerace, D. Valente, S. Montangero, A. Auffèves, and M. F. Santos, EPL (Europhysics Letters) 106, 54003 (2014).
  * Marcos-Vicioso  _et al._ (2018) A. Marcos-Vicioso, C. López-Jurado, M. Ruiz-Garcia, and R. Sánchez, Phys. Rev. B 98, 035414 (2018).
  * Bhandari  _et al._ (2021) B. Bhandari, P. A. Erdman, R. Fazio, E. Paladino, and F. Taddei, Phys. Rev. B 103, 155434 (2021).
  * Iorio  _et al._ (2021) A. Iorio, E. Strambini, G. Haack, M. Campisi, and F. Giazotto, Phys. Rev. Applied 15, 054050 (2021).
  * Balachandran  _et al._ (2019) V. Balachandran, S. R. Clark, J. Goold, and D. Poletti, Phys. Rev. Lett. 123, 020603 (2019).
  * Kargı  _et al._ (2019) C. Kargı, M. T. Naseem, T. Opatrný, O. E. Müstecaplıoğlu, and G. Kurizki, Phys. Rev. E 99, 042121 (2019).
  * Ordonez-Miranda  _et al._ (2017) J. Ordonez-Miranda, Y. Ezzahri, and K. Joulain, Phys. Rev. E 95, 022128 (2017).
  * Aligia  _et al._ (2020) A. A. Aligia, D. P. Daroca, L. Arrachea, and P. Roura-Bas, Phys. Rev. B 101, 075417 (2020).
  * Upadhyay  _et al._ (2021) V. Upadhyay, M. T. Naseem, R. Marathe, and O. E. Müstecaplıoğlu, Phys. Rev. E 104, 054137 (2021).
  * Díaz and Sánchez (2021) I. Díaz and R. Sánchez, New Journal of Physics 23, 125006 (2021).
  * Tesser  _et al._ (2022) L. Tesser, B. Bhandari, P. A. Erdman, E. Paladino, R. Fazio, and F. Taddei, New Journal of Physics 24, 035001 (2022).
  * Zhang  _et al._ (2017) Y. Zhang, X. Zhang, Z. Ye, G. Lin, and J. Chen, Applied Physics Letters 110, 153501 (2017).
  * Zhang  _et al._ (2018) Y. Zhang, Z. Yang, X. Zhang, B. Lin, G. Lin, and J. Chen, EPL (Europhysics Letters) 122, 17002 (2018).
  * Ghosh  _et al._ (2012) A. Ghosh, S. S. Sinha, and D. S. Ray, Phys. Rev. E 86, 011138 (2012).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * Carmichael (2002) H. J. Carmichael,  _Statistical Methods in Quantum Optics 1_ (Springer-Verlag, Berlin Heidelberg, 2002).
  * Joulain  _et al._ (2016) K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Katz and Kosloff (2016) G. Katz and R. Kosloff, Entropy 18 (2016).
  * Jiang  _et al._ (2015) J.-H. Jiang, M. Kulkarni, D. Segal, and Y. Imry, Phys. Rev. B 92, 045309 (2015).
  * Goury and Sánchez (2019) D. Goury and R. Sánchez, Applied Physics Letters 115, 092601 (2019).
  * Liu  _et al._ (2022) Y.-Q. Liu, D.-H. Yu, and C.-S. Yu, Entropy 24 (2022).
  * Gupt  _et al._ (2022) N. Gupt, S. Bhattacharyya, B. Das, S. Datta, V. Mukherjee, and A. Ghosh, Phys. Rev. E 106, 024110 (2022).
  * Gelbwaser-Klimovsky  _et al._ (2015) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki,  Advances In Atomic, Molecular, and Optical Physics, 64, 329 (2015).
  * Spohn (1978) H. Spohn, Journal of Mathematical Physics 19, 1227 (1978).
  * Kosloff (2013) R. Kosloff, Entropy 15, 2100 (2013).
  * Deffner and Campbell (2019) S. Deffner and S. Campbell,  _Quantum Thermodynamics_, 2053-2571 (Morgan and Claypool Publishers, 2019).
  * Ghosh  _et al._ (2019) A. Ghosh, V. Mukherjee, W. Niedenzu, and G. Kurizki, The European Physical Journal Special Topics 227, 2043 (2019).
  * Gupt  _et al._ (2021) N. Gupt, S. Bhattacharyya, and A. Ghosh, Phys. Rev. E 104, 054130 (2021).
  * Sinha  _et al._ (2013) S. S. Sinha, A. Ghosh, and D. S. Ray, Phys. Rev. E 87, 042112 (2013).
  * Thierschmann  _et al._ (2015a) H. Thierschmann, R. Sánchez, B. Sothmann, F. Arnold, C. Heyn, W. Hansen, H. Buhmann, and L. W. Molenkamp, Nature Nanotechnology 10, 854 (2015a).
  * Hartmann  _et al._ (2015) F. Hartmann, P. Pfeffer, S. Höfling, M. Kamp, and L. Worschech, Phys. Rev. Lett. 114, 146805 (2015).
  * Thierschmann  _et al._ (2015b) H. Thierschmann, F. Arnold, M. Mittermüller, L. Maier, C. Heyn, W. Hansen, H. Buhmann, and L. W. Molenkamp, New Journal of Physics 17, 113003 (2015b).
  * Pfeffer  _et al._ (2015) P. Pfeffer, F. Hartmann, S. Höfling, M. Kamp, and L. Worschech, Phys. Rev. Applied 4, 014011 (2015).
  * Schaller  _et al._ (2009) G. Schaller, G. Kießlich, and T. Brandes, Phys. Rev. B 80, 245107 (2009).
  * Zedler  _et al._ (2009) P. Zedler, G. Schaller, G. Kiesslich, C. Emary, and T. Brandes, Phys. Rev. B 80, 045309 (2009).
  * Wijesekara  _et al._ (2021) R. T. Wijesekara, S. D. Gunapala, and M. Premaratne, Phys. Rev. B 104, 045405 (2021).