Universal Behavior of the Coulomb-Coupled Fermionic Thermal Diode

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”.


Introduction
Heat management devices have attracted recent interest in nanoscale systems; they prevent overheating due to heat flow in desired areas of electronic circuits [1][2][3][4][5][6][7][8][9][10][11][12]. Following the theoretical proposal of a quantum dot thermal diode by Roukola and Ojanen [13], a number of studies have been carried out to achieve the rectification effect using simple quantum systems [14][15][16][17][18][19][20][21][22][23][23][24][25][26][27][28][29][30]. Yet, most of the works conducted to date rely on either temperature gradients of bosonic reservoirs or different coupling strengths between the system and the bath to break the inversion symmetry of the overall system [30]. For example, Werlang et al. [19] explored heat transport under the influence of strong coupling between two spins interacting with their respective bosonic baths. Miranda et al. [26] identified similar diode characteristics in the presence of different excitation frequencies between the coupled spins. Based on two-and four-terminal quantum dot setups, Tesser et al. [30] recently explored the roles of level degeneracy and temperature bias. However, much less attention is paid to achieving the rectification effect by means of the 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 independent of the system energy levels, interaction strength, bath spectrum, its temperature, and chemical potentials.
To showcase our findings, we considered a Coulomb-coupled quantum dot system and studied the interaction with two fermionic reservoirs with different temperatures and chemical potential. We find that, while the temperature gradient governs the overall heat flow direction, the magnitude of the heat current is primarily controlled by the chemical potential gradient. In contrast to the bosonic counterpart, the fermionic rectifier allows us to control the heat current and switching effects in a much more efficient way, even in the presence of tiny temperature differences. Most remarkably, we identified the universal nature of complete, partial, and no rectification conditions that are valid for all Coulombcoupled fermionic diodes.
The present work is organized as follows: We introduce the model and dynamics in Section 2, and the steady-state heat current in Section 3. The microscopic picture behind the thermodynamically consistent heat flow direction is summarized in Section 4. Universal characteristics based on two dimensionless parameters are presented in Section 5, and the role of an efficient thermal switch and the ideal rectification effects are discussed in Section 6. Finally, we conclude in Section 7.

Model and Dynamics
Our model consists of two quantum dots (QDs) that 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 References [13,27,[30][31][32], who recently introduced the study of thermal diodes and transistor effects in a wide variety of Coulomb blockade quantum dot devices. We further assume that each QD is tunnelcoupled to its fermionic reservoir [33]. While the electron transports between QDs are forbidden due to the Coulomb blockade [31], electron tunneling between QDs and their respective reservoirs permit heat flow from one reservoir to another through the coupled QD system. Within sequential tunneling [32] under the Coulomb blockade regime, each QD can only have two levels with the occupation number as either zero or one. The two QDs, as well as the temperature and chemical potential of the baths to which they are connected, are labeled by indices L and R [ Figure 1]. The Hamiltonian of the coupled QDs is then given by, Here, ε α is the lowest single-particle energy of QDs. Without loss of any generality, we further assume ε L < ε R . Here, U αβ is the positive Coulomb interaction energy between the electrons in different QDs, and d † α (d α ) denotes the creation (annihilation) operator for the α-th QD, whose eigenstates are |0 and |1 with eigenvalues 0 and ε α , respectively. Since the interaction energy in Equation (1) is diagonal in the eigenbases of the individual QDs, the 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 = |00 , |2 = |10 , |3 = |01 , |4 = |11 . For details, please refer to References [13,31,32]. The Hamiltonian of the α-th fermionic reservoir is defined as H α [30], where ε k is the energy of the non-interacting reservoir electron, the continuous wavenumber k, µ α is the chemical potential, and c † (c) represents the creation (annihilation) operator of the electron reservoir. The coupling between QDs and the respective reservoir is described by the tunneling Hamiltonian, where t αk is the tunneling amplitude. Interaction within sequential tunneling approximation imposes restrictions on simultaneous tunneling of more than one electron at a time [13,27,30,31]. Consequently, there are, in total, four authorized transitions: the left reservoir (L) induces transitions between 1 ↔ 2 and 3 ↔ 4, while the right bath (R) drives transitions between 1 ↔ 3 and 2 ↔ 4. We define transition energies ω ij = i − j , for i > j, where k is the eigenvalue of H D for the eigenstate |k . In the present case, they read as ω 21 = ε L , ω 42 = ε R + U, ω 43 = ε L + U, and ω 31 = ε R . The rates at which the above transitions occur are computed using the Lindblad master equation [34,35]. We implement the strong-coupling formalism following References [19,[36][37][38][39][40][41] to arrive at the master equation describing the time evolution of the density matrix ρ of the coupled QDs system (See Appendix A) under Born, Markov, and the secular approximation [19,21,34]. It is important to note here that the strong coupling formalism refers to the coupling between the dots, while the systembath coupling is still assumed to be weak so that the Born-Markov approximation can safely be implemented [19,[36][37][38][39][40][41]. This implies that the Lindbladians are obtained on the basis of the eigenstates of the full system Hamiltonian H D . Thus, the 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.

Evaluation of Steady State Heat Current
In the present model, particles cannot be exchanged between QDs as they interact only through the long-range Coulomb force. As a result, there is no particle flow in between reservoirs through the QD system [13,27,[30][31][32]. So, energy is exchanged only in the form of heat [See Figures 2 and 3]. The expression of the heat current can then be obtained using the standard procedure [42] starting from the von Neumann entropy of the system, defined as S[ρ(t)] = −k B Tr[ρ(t) ln ρ(t)]. Upon taking the time derivative of the von Neumann entropy, one obtains where we have used the master Equation (2) and Tr[ρ] = 0. The Sphon inequality [43] in the form of the second law of thermodynamics [44] for any Lindblad superoperator L can be written as Tr{L[ρ(t)](ln ρ(t) − ln ρ ss )} 0. Here, ρ ss is the steady-state population of the system, satisfying L[ρ ss ] = 0. Denoting the stationary state of L α , as ρ α ss , above, the inequality can be applied to both terms of the sum in Equation (3), yielding From Equations (3) and (4), we can write The above equation can be compared with the dynamical version of the second law given by [42,44] d dt This allows us to identify the heat current (energy flow rate) J α Q (t) associated with the α-th reservoir as J α For the steady state operation, the master Equation (2) drives the system towards a Gibbs-like stationary state, characterized by [41,42] For the detailed derivation of the steady state's heat current, we refer readers to the review articles by Kosloff and Gelbwaser et al. [42,44]. Inserting Equation (8) into Equation (7), we obtain the general expression for the heat current, which, under a steady state condition, simplifies to Here, the net decaying rate Γ α ij from |i to |j (i > j) is denoted as The first term represents the emission and the second term corresponds to absorption; γ α is the bare tunneling rate between the dots and respective reservoirs [ Figure 1] corresponding to the transition energy ω ij = i − j between eigenstates |i and |j controlled by the α-th reservoir. The tunneling of electrons into or out of QDs is primarily governed by FDF. In our model, there are four allowed transitions and both lead guides to two transitions each: the L lead drives transitions between |1 ↔|2 and |4 ↔|3 , while the R lead controls |1 ↔|3 and |4 ↔|2 transitions. For the sake of convenience of our analysis, the corresponding FDFs are expressed in terms of two dimensionless parameters, the effective tunneling barrier (χ α = (ε α − µ α )/U) and dimensionless thermal energy (ξ α = k B T α /U), as follows Now, our task is to first evaluate the full expression of Γ α ij at the steady state, and then find out the expression of the heat current. Under the steady state condition, the master Equation (2) is characterized byρ ss = 0, which reduces tȯ Thus, at the steady state, all net transition rates become equal to Γ, i.e., The four sets of equations in Equation (13) are not independent since ∑ i ρ ii = 1, which uniquely solves all of the state occupation probabilities as well as the heat current in terms of a single quantity Γ: (9), we can then evaluate the general expression of the heat current at the steady state as follows Thus, we finally arrive at the explicit analytical expression of the steady state heat current as To find out Γ, we rewrite Equation (13) in terms of f 1 L(R) and f 2 defined through Equation (12) and find out the steady state populations subject to the condition Using Equation (16) and Equation (17), we can construct where, Solving the above equation, the expressions for the steady state populations ρ ii are where |M| stands for the determinant of the matrix M. Using the above expressions, we can evaluate the final expression of Γ as follows Using Equation (20), the final expression of the steady state heat current in terms of dimensionless parameters χ α and ξ α can be written as where, This is the exact analytical expression of heat current derived under the Born-Markov master equation. Since U, γ L , γ R and X are all positive, it is immediately clear from Equations (21) and (22) that if ξ R > ξ L (T R > T L ), heat will flow from right to left (J R > 0) in accordance with Equation (15).

Microscopic Description of Heat Flow
Classically, heat flows according to the laws of thermodynamics, i.e., from high to low temperatures. Quantum mechanically, the flow of the heat current must be governed by a microscopic description, without violating the ultimate laws of thermodynamics [44][45][46][47][48]. Since the QDs are strongly coupled with each other, they form a four-level system as depicted in Figure 2, where the inter-dot transition is restricted due to the long-range Coulomb force. So, the energy transfer is only caused by the heat baths via the coupled states of the overall system [13,27,[30][31][32]. For T R > T L , the natural heat flow direction will be from R to L with J R Q > 0. In this case, the first excitation from the ground state must be guided by the cold bath, instead of a hot bath, which may appear paradoxical at first sight but makes the net transition rate Γ > 0 as required by Equation (15). If the cycle follows the opposite path 1 → 3 → 4 → 2 → 1 [ Figure 2(aI)], then J L Q > 0, as the UΓ 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 thermodynamics, the transition cycle must run in 1 → 2 → 4 → 3 → 1, initiated by the cold bath. This is seemingly paradoxical in the sense that, classically, we expect that, during the heat flow, energy is supplied by the hot bath and dumped into the cold bath. For the specific example, it is more favorable for the cold bath to make the 1 → 2 transition in Figure 2(aII), which costs ε L amount of energy than the 3 → 4 transition in Figure 2(aI), which requires (ε L + U) amount of energy. It is interesting to note that, for T L > T R , (i.e., J L Q > 0, the heat flows from the left to the right following Figure 2(bII); the first transition is still mediated by the cold bath between 1 → 3, as it requires less energy (ε R ) than the 2 → 4 transition (ε R + U) in Figure 2(bI).
It is important to emphasize that although we have taken ε L < ε R to draw the schematic energy level diagram of Figure 2, the magnitude of the heat current, in particular, depends only on the two dimensionless parameters {χ α , ξ α } via Equation (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, Figure 3 shows that the basic principle behind the thermodynamically consistent transition cycle remains the same even in the case of ε L = ε R .

Universal Characteristic Due to Magic Means
While the temperature gradient dictates the overall heat flow direction, chemical potential plays a very important role in determining the heat current's magnitudes. To illustrate, we note that FDFs f 1 L(R) and f 2 L(R) in Equation (12), expressed in terms of dimensionless thermal energy ξ L(R) = k B T L(R) /U and the effective tunneling barrier χ L(R) = (ε L(R) − µ α )/U, are constrained by 0 ≤ f 2 L(R) < f 1 L(R) ≤ 1 and become 0.5 at χ L(R) = 0 and −1, which correspond to µ L(R) = ω 21(31) = ε L(R) and µ L(R) = ω 43(42) = ε L(R) + U, respectively. According to Equation (10), if 0 ≤ { f 1 L(R) , f 2 L(R) } ≤ 0.5, it favors de-excitation and, equivalently, excitation is favored when 0.5 ≤ { f 1 L(R) , f 2 L(R) } ≤ 1, for the corresponding transition. This is unique to fermionic reservoirs and allows us to implement a systematic analytical scheme solely based on the FDFs, and is in sharp contrast to its bosonic counterpart. The span of f 1 L(R) and f 2 L(R) as function χ L(R) , can thus be classified into five domains D [1−5] [ Figure 4: Main] and the total FDF, defined as f L(R) = f 1 L(R) + f 2 L(R) , is found to be spread between 0 ≤ f L(R) ≤ 2 [ Figure 4: Inset].  ; span of D 1 [5] and D 2 [4] are not fixed and strongly depend on the value of ξ L(R) . The color gradients signify the magnitude of |J Q |.
hence, excitation and de-excitation are favored for all four transitions yielding non-zero heat currents. On the contrary, domain D 3 , parameterized by −1 ≤ χ L(R) ≤ 0, is sharply defined between the transition points f 2 L(R) and f 1 L(R) . As opposed to D 3 , both f 1 L(R) and f 2 L(R) are closer to 1[0] in D 2 [4] , so that the heat flux decreases in D 2 [4] relative to D 3 . In D 3 , both f 1 L(R) and f 2 L(R) take values, such that excitation from |1 → |2 (|2 → |4 ) and de-excitation from |4 → |3 (|3 → |1 ) can occur simultaneously at optimal rates. The condition becomes ideal at the midpoint of D 3 (point B in Figure 4), which corresponds to χ L(R) = −0.5. To find out, analytically, the optimal value of χ L(R) for which |J Q | becomes maximum, first we have to differentiate J Q ≡ J R Q (Cf. Equation (21)) w.r.t. χ L(R) and set that equal to zero. Now, by differentiating Equation (21) w.r.t χ R , for fixed {χ L , ξ L , ξ R }, we obtain, As, U > 0 and γ L , γ R = 0, Equation (23) implies that the only a non-trivial solution for maximizing |J Q | is equivalent to maximizing Γ, which is given by the criteria It is clear from Equation (23) that the heat current vanishes under two limiting conditions: Again, from Equation (12), we can write So, Equation (25) signifies that, if f 1 is certainly 1. In both cases, the heat current vanishes. Thus, the only condition for the nonzero heat current reduces to f 1 R = 0 and f 2 R = 1. Under these conditions, after solving Equation (24), we obtain the value of χ R = −0.5, for which |J Q | attains the maximum: So, the heat current will be maximum when χ R = −0.5 [ Figure 5a]; similarly, when χ L varies, χ L = −0.5 is the criteria for having the highest heat current. Hence, we can conclude that |J Q | is maximum when both χ L(R) = −0.5, irrespective of ξ L(R) , which is supported by numerical simulations [ Figure 5b,c] and also follows from analytically derived conditions f L(R) = 1 [Cf. Equation (24)] for the maximum heat current [ Figure 4: Inset]. Now, putting these conditions in Equation (21), we can evaluate the exact analytical expression of the maximum heat current as Several remarks are now in order: • The maximum heat current is obtained at χ L(R) = −0.5 and the magnitude only depends on the temperature of two leads, their tunneling rates, and the Coulomb interaction between the dots. Since −0.5 is the mean of 0 and −1, which are transition points of f 1 L(R) and f 2 L(R) , respectively (point B in Figure 4: main), and this mean is also independent of other controlling parameters, we term this number "−0.5" as the "universal magic mean". The significance of the point B lies in the fact that at the magic mean χ L(R) = −0.5, the chemical potential of the left (right) lead becomes exactly equal to the means of the transition energies ω 21(31) and ω 43(42) driven by that bath: and, therefore, provides maximum control over the L(R) bath to guide both the absorption and decay simultaneously at the maximum Γ, with the resulting maximum |J Q |. • With the increase of ξ L(R) , |J Q | also spreads out, keeping the maxima point fixed at χ L(R) = −0.5. This precisely indicates that the maximum heat current will invariably be obtained at χ L(R) = −0.5, irrespective of all ξ L,R [ Figure 5b,c]. The 3D plots along χ L are squeezed for smaller values of ξ L(R) [ Figure 5b] (assuming ξ L < ξ R , as in Figure 2), and if we increase ξ L (ξ R ), it expands along both positive and negative χ L(R) , leaving out the point of maxima intact at χ L(R) = −0.5 [ Figure 5c]. This fundamental feature of the "magic mean" makes it truly universal.

Efficient Heat Current Modulator
With the increase of ξ L(R) , D 1 [5] becomes narrower while D 2 [4] spreads out, without affecting D 3 [ Figure 4]. As a consequence, |J Q | spreads out with ξ L(R) , keeping the point of maxima fixed at χ L(R) = −0.5. Now, if ξ L(R) 1, then the span of D 2 [4] is very small compared to D 1 [5] , while domain D 3 always remains in between −1 ≤ χ L(R) ≤ 0. So, there are effectively two domains: (I) Domain ON, where |J Q | = 0 and (II) Domain OFF, where |J Q | = 0 and the switching effect becomes more prominent if ξ L(R) 1. Thus, we can operate our model as an efficient thermal switch to the on-off heat current just by shifting the domain from ON to OFF through the change of χ L(R) or, in turn, the controllable experimental parameter ε L(R) [ Figure 6a]. This is the basic underlying principle behind the switching effect for each Coulomb-coupled fermionic thermal diode. It becomes more prominent for a smaller cold bath temperature ξ = k B T/U and can be achieved more easily by varying U, instead of lowering T to a smaller value.
Finally, the model operates as an efficient thermal diode with the rectification factor (R) approaching 1, even in the presence of a low-temperature gradient. This produces a clear advantage over all previously proposed models with the R factor defined as [4,16,18,22,36], where ∆ξ = ξ R − ξ L is identified as the temperature gradient, for ξ R > ξ L . With ∆ξ > 0, i.e., ξ R(L) ≡ ξ hot(cold) , the heat current flows from R to L, fulfilling J R Q > 0, and if we exchange temperatures of the bath, then the temperature gradient becomes −∆ξ or ξ R(L) ≡ ξ cold(hot) and, consequently, heat flows from L to R, satisfying J L Q > 0. Now, if the heat current vanishes upon reversing the temperature gradient, i.e., |J Q | is finite in one direction but null in the other, then complete rectification is achieved with R → 1, irrespective of the value of ∆ξ, whereas R → 0 corresponds to no rectification, i.e., no change in |J Q | upon inverting the temperature gradient. So, the current asymmetries in two directions are the primary criteria behind positive rectification: (i) It is clear from Figure 4 that for a given ∆ξ, |J Q | depends on |∆χ|. If χ hot = χ cold , or ∆χ = χ hot − χ cold = 0, the heat currents will be symmetric in both ways; therefore, R = 0. (ii) From Figure 6a, we find that the variation of |J Q | w.r.t χ α , is completely symmetric about the magic mean −0.5, i.e., |J Q (χ α )|=|J Q (−1 − χ α )|. As a result, following Equation (29) and Figure 6a, we find if χ hot = −1 − χ cold orχ = 1 2 (χ hot + χ cold ) = −0.5, then J R Q (∆ξ) = J L Q (−∆ξ), yielding R = 0. Thus, R can be zero; either (i) the difference between the effective tunneling barriers ∆χ = 0, or (ii) the mean effective tunneling barrier is equal to the magic mean (χ = −0.5); otherwise, rectification occurs.
Once |∆χ| = 0 andχ = −0.5, effective tunneling barriers χ L(R) make dissimilar effects on |J Q | upon reversing the temperature gradient and, consequently, R becomes nonzero. For a fixed value of |∆ξ|, R increases as ∆χ shifts from zero and approaches one asχ significantly deviates from the magic mean [ Figure 6b]. The reason behind this, for smaller |∆χ|, the effect of tunneling barriers on |J Q | are comparable upon reversing the temperature gradient, yielding partial rectification (0 < R < 1). However, with the increase in |∆χ|, the effect of χ L(R) is no longer compatible and it creates larger heat current asymmetry between the two directions, leading to complete rectification behavior with R → 1. As the heat current is symmetric aboutχ = −0.5, we only consider positive variations ofχ in Figure 6b. Moreover, variations of R with |∆χ| depend on the value of |∆ξ|; larger ∆ξ implies better rectification. As ξ cold is less than ξ hot and the transition cycle is always initiated by the cold bath, the heat current executes stronger dependence on ξ cold than ξ hot for a given change in |∆ξ|. Thus, the variation of R with |∆χ| yields an appreciable change when a given |∆ξ| is altered due to the change in ξ cold than ξ hot [ Figure 6c: Main]. So, the complete rectification is more favorable when ξ cold 1 [ Figure 6a]. We may have complete rectification even without considering the absolute temperature of the cold bath close to zero, as we can vary ξ cold = k B T/U by changing U and keeping T finite. Finally, Figure 6c (Inset) shows that heat currents become ∼ 10 −6 times smaller upon reversing the temperature gradient, which corresponds to the R ≈ 1 curve in Figure 6b.

Conclusions
To conclude, the present model can be implemented as an efficient thermal switch as well as a thermal rectifier to modulate the 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 the maximum heat current at the magic mean "−0.5" in terms of dimensionless physical parameters, is the major finding 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 parameter, and it truly reflects the impact of chemical potential to decide the magnitudes of the heat current. The present protocol is unique to fermionic systems and can be applied to more complicated three or multi-terminal fermionic devices. The straightforward generalization of the present scheme to the multi-terminal set-up would involve 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 interactions, and can be simultaneously tunnel-coupled to their respective reservoirs. Their discrete energy levels provide energy-selective transport and can be tuned via the application of the external gate voltages. In view of the recent experimental advances in Coulomb-coupled quantum-dot systems [49][50][51][52], our findings will have important implications in designing novel thermal devices and opening up potential applications in controlling the thermal current at the nanoscale level.