Heat rectification through single and coupled quantum dots

We study heat rectification through quantum dots in the Coulomb blockade regime using a master equation approach. We consider both cases of two-terminal and four-terminal devices. In the two-terminal configuration, we analyze the case of a single quantum dot with either a doubly-degenerate level or two non-degenerate levels. In the sequential tunneling regime we analyze the behaviour of heat currents and rectification as functions of the position of the energy levels and of the temperature bias. In particular, we derive an upper bound for rectification in the closed-circuit setup with the doubly-degenerate level. We also prove the absence of a bound for the case of two non-degenerate levels and identify the ideal system parameters to achieve nearly perfect rectification. The second part of the paper deals with the effect of second-order cotunneling contributions, including both elastic and inelastic processes. In all cases we find that there exists ranges of values of parameters (such as the levels' position) where rectification is enhanced by cotunneling. In particular, in the doubly-degenerate level case we find that cotunneling corrections can enhance rectification when they reduce the magnitude of the heat currents. For the four-terminal configuration, we analyze the non-local situation of two Coulomb-coupled quantum dots, each connected to two terminals: the temperature bias is applied to the two terminals connected to one quantum dot, while the heat currents of interest are the ones flowing in the other quantum dot. Remarkably, in this situation we find that non-local rectification can be perfect as a consequence of the fact that the heat currents vanish for properly tuned parameters.


Introduction
Rectification is the phenomenon for which the magnitude of a current flowing in a system depends on the sign of the bias applied. In other words, by reversing the bias the current not only changes sign but also its magnitude. Perfect rectification is obtained when the current can flow only in one direction. The most familiar example of rectification is the one occurring in diodes, a two-terminal electronic component which allows the flow of charge current primarily in one direction, i.e. presenting low resistance in one direction and high resistance in the other.
Thermal rectification, i.e. the rectification of heat currents, occurs in a two-terminal system when the absolute value of the heat flux changes by reversing the sign of temperature bias applied to the two leads. This phenomenon has recently attracted increasing interest as a mean to improve thermal management in nanoscale systems, for example by blocking the flow of heat in certain areas of an electronic circuit to prevent overheating. Such interest is fueled by recent advancements in the experimental realization of nanostructured devices where thermal fluxes can be measured [1][2][3][4][5][6]. Thermal rectification was first observed experimentally a long time ago in reference [7]. More recently, in solid-state quantum systems it has been theoretically studied in references . In electronic nanoscale systems thermal rectification has been studied theoretically in hybrid quantum devices [15,16,22,31,32,34] and quantum dots (QDs) [18,20,30,[36][37][38][39], and experimentally measured in references [40][41][42][43][44]. In reference [33] thermal rectification has been calculated for a multi-level bosonic quantum system consisting of a nonlinear resonator attached to two baths.
Focus of the present work is to investigate how heat can be rectified using QD-based devices. As mentioned above, few papers on heat rectification in QDs are available in the literature. In reference [36], motivated by the experiment reported in reference [41], thermal rectification has been investigated for a two-level QD using a nonequilibrium Green function method and focusing on the role of the energy-dependence of the tunnel couplings between QD and leads. In references [18,20,39] the case of multiple capacitively-coupled QDs was considered: while in reference [18] all QDs were connected to two leads, in references [20,39] each of the two QDs considered were attached to one lead only, so that heat can be transported only by electronic fluctuations. In reference [37] a mean-field approximation was used to calculate self-consistently the heat rectification of a single QD. The role of interference, quantum superposition and level degeneracy on the heat rectification was studied in reference [30] for various systems of coupled QDs, using the master equation approach up to sequential tunneling processes. Finally, the case of an array of QDs was studied in reference [38] using the Keldysh Green's function technique.
An important ingredient to rectify heat is the presence of non-linearities in the spectrum (e.g. due to the combined effect of confinement and electron-electron interaction). This fact can be easily understood by noticing that, in the absence of interactions, the heat current can be calculated using the Landauer-Büttiker approach, i.e. by an energy integral of the transmission probability of the QD multiplied by the difference of the Fermi distribution functions of the two terminals. Since the temperatures enter only the distribution functions (indeed the transmission probability consists of a set of narrow Lorentzian functions of energy, one for each QD level), an inversion of the temperature bias simply gives rise to a change of sign of the current, thus no rectification. Non-linearities, however, are not enough. Indeed, to obtain rectification a necessary condition is to break the mirror-symmetry of the system, for example by coupling the system to the left and to the right terminals by a different extent.
In this paper we study heat rectification for QDs in the Coulomb blockade regime using the master equation approach [20,30], with a particular emphasis on the role of second-order cotunneling contributions. We consider both the two-terminal and the four-terminal configurations (see figures 1(a)-(b) and (c)-(d), respectively). For the former, we analyze the case of a single QD with either a doubly-degenerate level or two non-degenerate levels. The two reservoirs, labelled with L and R in figures 1(a) and (b), are characterized by their temperatures (T L = T + ΔT/2 and T R = T − ΔT/2) and their chemical potentials (μ L = Δμ/2 and μ R = −Δμ/2). We are interested in studying the heat current that flows through the system when a temperature bias is applied between the reservoirs. Furthermore, we assume that no work is performed on the system and consider both the open-circuit setup (where the charge current vanishes, see figure 1(e)) and the closed-circuit setup (where the bias voltage is set to zero, see figure 1(f)). For the sake of definiteness, we focus on the heat current flowing between the left terminal and the QD. Referring to figures 1(a) and (b), we define the forward heat current J + as the one relative to ΔT > 0 and the backward heat current J − the one relative to ΔT < 0. The laws of thermodynamics assure that heat will flow from left to right if ΔT > 0 (forward bias), or from right to left if ΔT < 0 (backward bias). When |J + |, induced by a forward bias, is different from |J − |, induced by a backward bias, we define the heat rectification coefficient as The definition is such that |R| 1. In particular, R = 0 means that no rectification takes place, while R = ±1 means that we have perfect rectification (i.e. the heat current is finite in one direction, and null in the other). We first analyse the sequential tunneling regime (assuming up to single occupancy of the QD) deriving, when possible, analytical expressions for the heat currents. Remarkably, we could derive an upper bound for rectification in the closed-circuit setup with a doubly-degenerate level, and prove the absence of a bound for the case of two non-degenerate levels. We also analyze the behaviour of currents and rectification as functions of the levels' position and the temperature bias. The most important part of the paper deals with the effect of cotunneling contributions, including both elastic and inelastic processes, in all setups of the two-terminal configuration. The most remarkable results are the following: (i) in the open-circuit setup of the doubly-degenerate level case, cotunneling yields a finite (though little) rectification, contrary to what Figure 1. Sketch of the systems considered: grey circles represent QDs, while red, blue and green objects represent the reservoirs. Panels (a) and (b) refer to the two-terminal setup, while panels (c) and (d) refer to the four-terminal setup. For the two-terminal setup, the left (right) reservoir is characterized by a temperature T L = T + ΔT/2 (T R = T − ΔT/2) and a chemical potential μ L = Δμ/2 (μ R = −Δμ/2). Panel (a) represents the forward bias configuration, where ΔT > 0 (the left reservoir is hot and the right reservoir is cold), while panel (b) represents the backward bias configuration, where ΔT < 0 (the left reservoir is cold and the right reservoir is hot). J ± is the heat current flowing in the left lead in the forward (backward) bias configuration. In the four-terminal setup, the two QDs (identified by the arrows ↑ and ↓) are each connected to two reservoirs and are Coulomb-coupled to each other (with charging energy E C ). We refer to L1, L2 and QD ↑ as the drive circuit, while to R1, R2 and QD ↓ to the drag circuit and we set T L1 = T + ΔT/2, T L2 = T − ΔT/2, T R1 = T and T R2 = T. All reservoirs are kept at the same chemical potential. Panel (c) depicts the forward bias configuration, where ΔT > 0, while panel (d) depicts the backward bias configuration, where ΔT < 0. We are interested in the heat currents J ± R1 and J ± R2 flowing in the lead R1 and R2, respectively. Panels (e) and (f) illustrate the open-circuit and the closed-circuit setups, respectively. In the former case no charge current flows through the system, while in the latter the two reservoirs are electrically connected (μ L = μ R ). happens when only sequential processes are considered; (ii) in the closed-circuit setup of a doubly-degenerate level, cotunneling corrections to the forward heat current are opposite to the corrections to the backward heat current, thus yielding rectification enhancement when cotunneling lowers the magnitude of the two currents (an analogous result was reported in reference [33]); (iii) in the case of two non-degenerate levels, in the open-circuit setup, cotunneling always increases the currents with respect to the sequential regime; (iv) in all cases there exists ranges of values of the levels' position where rectification is enhanced by cotunneling.
For the four-terminal configuration, we analyze the case of two Coulomb-coupled QDs, each connected to two terminals (see figures 1(c) and (d)). Such a setup has been actually realized in references [45][46][47][48][49][50][51][52][53][54][55]. This is a non-local configuration, where the temperature bias is applied to the terminals L1 and L2 on the left (drive circuit), while the heat currents of interest are the ones flowing in terminals R1 and R2 on the right (drag circuit). Remarkably, in this situation we find that non-local rectification (defined for the currents in the drag circuit) can reach the ideal value as a consequence of the fact that the heat currents in the drag circuit can change sign (thus going to zero) as a function, for example, of the energy level of one of the QDs. The absolute value of the heat currents in the drag circuit, however, is small when compared to the heat currents in the drive circuit.
In addition, we consider the case where the tunnel couplings between QDs and leads depend on energy, since this situation occurs in experimental realizations [53]. We find that the heat currents in the drag circuit have similar amplitude but opposite signs, meaning that if heat is extracted from one reservoir, a similar amount of heat is deposited into the other.
The paper is organized as follows: in section 2 we define the systems under investigation and we describe the model Hamiltonian, while in section 3 we discuss the results we obtain. We first consider the sequential tunneling regime for a single QD with a doubly-degenerate level in section 3.1.1, with two non-degenerate levels in section 3.1.2. Then we discuss the results obtained when cotunneling contributions are accounted for in section 3.2 (with a doubly-degenerate level in section 3.2.1, and with two non-degenerate levels in section 3.2.2). Section 4 is devoted to the results obtained with the four-terminal (non-local) configuration, with two QDs each coupled to two reservoirs. The summary can be found in section 5. Details of the calculations relative to the master equation in the sequential tunneling regime are reported in appendix A, while cotunneling contributions are reported in appendices B and C.

System and model
We consider a system consisting of a QD with two levels (relative to spin up and spin down), whose Hamiltonian reads where ↑ = − Δ /2 and ↓ = + Δ /2 are the energy of the two levels. Heren σ =ĉ † σĉσ is a number operator, whileĉ † σ andĉ σ are creation and destruction fermionic operators, respectively, for an electron with spin σ in the QD. The first two terms describe the discrete levels of the QD, while the last one accounts for the Coulomb repulsion between the electrons within the QD (E C represents the charging energy). We assume the spacing Δ to be much smaller than E C , so that the electrostatic interaction plays a fundamental role in the transport properties of the system.
The QD is tunnel-coupled to two electronic reservoirs characterized by a well defined temperature T α and chemical potential μ α , whose Hamiltonians are given by whereb αkσ andb † αkσ are, respectively, the destruction and creation operators for electrons in lead α = L, R with energy αkσ , spin σ and momentum k. The coupling Hamiltonian reads where t α is the tunneling amplitude between the QD and lead α. The Hamiltonian H T is not symmetric in the coupling, namely t L = t R . Indeed, this is the condition needed to obtain a finite rectification.
In the final part of the paper we explicitly consider a four-terminal (non-local) configuration consisting of two single-level QDs, each tunnel-coupled to two reservoirs as sketched in figures 1(c) and (d), whose Hamiltonian is given by equation (2) where σ specifies the QD (↑/↓ for the QD on the left/right). In this case the coupling Hamiltonian reads while the Hamiltonian for the reservoirs is given by equation (3), with the index α now taking the following values: L1, L2, R1, R2. The state of the QD (or QDs) is specified through the probability P({σ}) of finding the QD in the electronic configuration described by the set of occupancies {n σ } of its levels (with n σ = 0, 1). Tunneling processes change the state of the QD thus modifying the occupancies of the levels from one configuration to another. In what follows we shall first consider the sequential tunneling regime which accounts for tunneling processes involving single-electron hopping through the tunnel barriers representing the coupling between the QD and the leads. Therefore, the transition rates of such processes are obtained from the Fermi golden rule up to the leading order in the coupling Hamiltonian H T (see appendix B). The tunneling constants Γ α characterizing the interaction between the QD and lead α are defined as where D α is the density of states of lead α = L, R at the Fermi energy. In part of our analysis, we will further assume that the charging energy E C is the largest energy scale. This allows us to neglect all electronic configurations in which the total number of electrons in the QD exceeds one. Therefore we can describe the state of the QD by specifying the probability of finding the QD in the state with zero electrons P 0 , with one electron in a level with spin up (P ↑ ), and with one electron in a level with spin down (P ↓ ). The master equation needed to determine such probabilities and the expressions of the currents are reported in appendix A.

Results
Let us first fix the notation: we denote by I c and I, respectively, the charge and energy current entering the QD from the left lead. The heat current J flowing through the left lead is thus expressed as where −e is the electronic charge. Note that both in the open-circuit setup, where I c = 0 (figure 1(e)), and in the closed-circuit setup, where μ L = μ R = 0 (figure 1(f)), we have that heat and energy current coincide (J = I).
Before discussing in details our results in the various regimes and configurations, we first show the rectification coefficient obtained for two representative cases. Namely, for a single QD (see figures 1(a) and (b)) with two non-degenerate levels in the closed-circuit setup and for a pair of QDs in the four-terminal setup (see figures 1(c) and (d)). In figure 2(a) the rectification coefficient is plotted for the former case as a function of the energy of the level . While R includes cotunneling contribution (solid red curve), R seq accounts for sequential tunneling processes only (green dashed curve). Figure 2(a) shows that R is typically not very large (for the parameters used here, the maximum value of R is of the order of 2%) and can take both negative and positive values depending on the position of the levels of the QD. Remarkably, we find that cotunneling contributions can increase the rectification in a wide range of values of . All details will be discussed in section 3.2.2.
In figure 2(b) the rectification coefficient for the four-terminal (non-local) configuration is plotted as a function of the level ↑ of the QD on the left-hand-side. We consider two values of temperature bias: ΔT = 0.1T (red solid curve) and ΔT = 0.3T (blue dashed-dotted curve). Remarkably, the blue curve spans the whole range of values of R ([−1, 1]), while the red curve takes values in the range [−1, 0.4]. Overall, rectification is large in quite wide ranges of values of ↑ . As we will see in more details in section 4, however, the heat currents are rather small when compared with the single QD setup.
In the following sections we will describe the results obtained within the sequential tunneling regime (section 3.1) and the results obtained accounting for the cotunneling contributions (section 3.2).

Sequential tunneling regime
In this section we will assume that the charging energy E C is so large that we can neglect all electronic configurations in which the total number of electrons in the QD exceeds one. This assumption, which will be lifted in section 3.2, allows us to obtain analytical results. The results shown in the following sections 3.1.1 and 3.1.2 are in agreement and largely extend the results presented in reference [30]. In particular, on the one hand, we will identify upper bounds for rectification and, on the other, we will discuss the relevant mechanisms allowing for the optimization of rectification both for degenerate and non-degenerate levels.

Degenerate level
Let us consider the degenerate case Δ = 0 in which the charge current can be written as while the heat current, in accordance with reference [30], takes the form thus showing that the heat current is proportional to the charge current. A direct consequence of this is that in the open-circuit setup (where there is no charge flow) the heat current is zero.
In the closed-circuit setup, however, the heat currents are finite and the rectification can be written as [30] showing that a necessary condition to obtain rectification is that Γ L = Γ R (this condition reflects the necessity to break the mirror-symmetry of the system). In figure 3 we plot the absolute values of the forward and backward heat currents and the resulting rectification. In particular, in panel (a) J + and |J − | are plotted as functions of the QD levels' energy . When is zero (i.e. when is aligned with the common chemical potential of the leads), both heat currents vanish because in this symmetric situation the sequential processes relative to the two leads cancel out. When | | 0, the heat currents decrease exponentially because the QD is locked in the same state and the electrons cannot tunnel. Indeed, P 0 goes quickly to 1 (and P ↑ = P ↓ goes to zero) when increases to positive values over the scale set by k B T, since electrons do not have enough energy to enter the QD, while P 0 goes quickly to zero (and the QD gets occupied) when decreases to negative values, as the energy level of the QD goes well below the chemical potential of the leads. Notice that the currents display two asymmetric maxima at | | ≈ 2.5k B T, and the maximum at ≈ −2.5k B T is lower than the one at ≈ 2.5k B T. The reason for this is that the probability of having one electron in the QD is higher when is negative (and of the order of k B T), as compared to when is positive, so that fewer electrons can enter the QD (thus contributing to the current) because the charging energy does not allow any other electron to tunnel in the QD.
In panel (c) the rectification is plotted as a function of the QD levels' energy . Like the currents, the rectification goes to zero when = 0 and when | | 0. Intuitively, we can understand that R gets suppressed for large positive values of by noting that in this situation Coulomb interaction plays little role (the QD is essentially unoccupied). Thus the QD virtually behaves as a non-interacting one where rectification does not occur. It turns out that the magnitude of rectification has two asymmetric maxima at | | ≈ 1.6k B T, where |R| ≈ 0.010 and |R| ≈ 0.015. It is worth stressing that both heat currents and rectification are close to their maximum when is within the interval 1.6-2.5k B T, so that the QD operates as a heat rectifier to the best of its capabilities. In panels (b) and (d), we plot the heat currents (J + and |J − |) and the rectification, respectively, as functions of the temperature bias ΔT. Increasing |ΔT|, the currents and their separation (J + + J − ) grow. Notice that when ΔT/T 1 and the currents are in the linear-response regime, the rectification vanishes (R = 0 at ΔT = 0). Increasing ΔT, the currents go beyond the linear-response regime and R varies linearly with ΔT. For values of ΔT/T larger than 0.5 the rectification is sublinear, but monotonous, thus reaching its maximum at |ΔT| = 2T.
Interestingly, in this configuration it is possible to find an upper bound that limits the rectification in this system. Let us consider the rectification parameter R as a function of f L and f R , see equation (10), and look for its maximum over the possible values that the Fermi distribution functions can take. It is important to notice that since f L and f R are evaluated at the same energy , the quantities f L ( ) and f R ( ) cannot take arbitrary values between 0 and 1. Actually, it is easy to see that both f L , f R 1/2, when 0, or both f L , f R 1/2, when 0 (recall that in the closed-circuit setup μ L = μ R = 0), regardless of the temperatures T L and T R . With this constraint taken into account, it is possible to prove that the maxima of the function |(f L − f R )/(2 + f L + f R )|, appearing in equation (10), occur when f L = 0 and f R = 1/2, or when f L = 1/2 and f R = 0. In particular, f L,R = 0 corresponds to /(k B T L,R ) 1, i.e. T L,R /k B , and f L,R = 1/2 corresponds to /(k B T L,R ) 0, i.e. T L,R /k B . By substituting the above values of f L and f R , we obtain the following upper bound

Non-degenerate levels
Let us now consider the non-degenerate case Δ = 0 in which the spin degeneracy of the level of the QD is broken, for example, through the Zeeman effect by applying a magnetic field. Closed-circuit setup. In this configuration we set μ L = μ R = 0. In figure 4 we plot the heat currents calculated as functions of the average levels' energies . We first note that the heat current (panel (a)) resembles the behavior found for the degenerate case, reported in figure 3, with important differences. Like in the degenerate case, section 3.1.1, the currents have two asymmetric maxima and are suppressed exponentially at large | | (the currents are suppressed when 0, since the QD is mostly empty, and when 0, the QD is mostly occupied by one electron). However, at = 0, the sequential tunneling processes from the left lead do not cancel out with the ones from the right lead because now the levels have different energies. As a consequence, when = 0 both heat currents are finite. Moreover, they present a local minimum when Δ /2, independently of the values of ΔT and tunneling constants. This can be understood by noting that, at least when Δ k B T L , k B T R , electron transport is mainly due to the lower level ( ↑ ), which is aligned with the common chemical potential of the leads, while the upper level is too high in energy, thus hardly populated (P ↓ 0). In this situation, however, energy current is minimum since where the levels are nearly degenerate. Notice that the distance between the two maxima is controlled by the average thermal energy k B T. The resulting rectification is plotted in panel (b) as a function of . For large values of | |, R behaves as in the degenerate case, while when | | < 5k B T the rectification oscillates between positive and negative values presenting an absolute maximum at ≈ Δ /4, close to the heat current minimum. Moreover, the negative dip on the left occurs approximately at −Δ /2, while the negative dip on the right occurs approximately at Δ : both correspond to values of heat currents between the local minimum and the maxima. The positions of such peaks and dips, however, depend also on the other parameters of the system. The behavior of heat currents and rectification with the bias ΔT is essentially the same as the one found in section 3.1.1. Notice that the upper bound found for the degenerate case, equation (11), does not apply here. For = k B T/2 and the parameters used in figures 3 and 4 we find R 0.11 for the largest value of ΔT, which is larger than the bound R max 0.067.
Intuitively one can expect the rectification to be optimized by maximizing the asymmetry between the tunneling constants Γ L and Γ R , and for large temperature bias. By using equation (A10) we can prove that in the closed-circuit setup it is possible to reach perfect rectification when the parameters satisfy the conditions: Such conditions are represented, for the forward configuration, in the energy diagram in figure 5 (left panel), where the lead L is hot and the lead R is cold. Because of (iii) and (iv), the distribution functions can be approximated as f L↑ ≈ f L↓ ≈ 1/2 and f R↑ ≈ 1 − f R↓ ≈ 1, so that the forward heat current takes the form (refer to appendix A for the notation). By taking the limit Γ L Γ R , one finds that Λ ≈ 3Γ 2 L /4 so that J + ≈ Γ R Δ /3. The physical origin of this expression for the heat current can be understood by looking at the energy diagram in figure 5 (left panel). Since the right lead's temperature is much smaller than | ↑ |, the tunneling rate of the process which transfer an electron from ↑ to lead R and the one which transfer an electron from lead R to ↓ are suppressed. Thus, the heat transport happens mainly through two following processes. The first one involves an electron tunneling from lead L to ↓ and, from there, to lead R. In the second one, the electron starts in the right lead, tunnels into ↑ , and then arrives in lead L. Such processes occur on the same typical time, of the order of 1/Γ R . Since they involve different QD levels, they transfer different amounts of heat.
8T. The right barrier is thicker than the left one because the tunneling constants satisfy Γ R Γ L . In the left panel the device is in the forward configuration, namely T L > T R , while, in the right panel the leads' temperatures are exchanged and the device is in the backward configuration. In both panels the chemical potentials of the leads (μ F ) are set to zero (dotted line). The rectification coefficient turns out to be R = 0.31.
For the backward configuration, where the lead L is cold and the lead R is hot, the energy diagram is represented in figure 5 (right panel). Because of (iii) and (iv), the distribution functions can be approximated as f L↑ ≈ 1 − f L↓ ≈ 1 and f R↑ ≈ f R↓ ≈ 1/2, so that the backward heat current takes the form By taking the limit Γ L Γ R , one finds that Λ ≈ Γ 2 L so that J − ≈ Γ R ↑ /2. Note that J − is, correctly, a negative quantity. Also in this case the physical origin of this expression for the heat current can be understood by looking at the energy diagram in figure 5 (right panel). The low temperature of the lead L suppresses the tunneling rate of the process in which an electron tunnels from ↑ to the lead L and of the one in which an electron tunnels from lead L to ↑ . Therefore, heat transport is dominated by the process in which an electron in lead L tunnels into ↑ and, from there, tunnels into lead R, thus transferring an amount of heat equal to ↑ in a typical time 1/Γ R . This gives rise to the current in equation (13). By plugging in the expressions for J + and J − into the definition of R, equation (1), we find R = (Δ /3 + ↑ /2)/(Δ /3 − ↑ /2). By imposing condition (ii) we find R 1. The drawback is that both heat currents are suppressed, since we have assumed a small value of Γ R .
Open-circuit setup. Although the charge current is zero, the fact that the two levels have different energy allows heat transport, contrary to what happens in the degenerate case (section 3.1.1). Indeed, to nullify the charge current (equation (A8)) the rate of electrons tunneling into the lower level (first term in the square bracket of equation (A8)) has to cancel out with the rate of electrons tunneling into the upper level (second term in the square bracket), namely In this condition the energy current is finite and reads where the last equality is obtained using equation (A5). In figure 6, the heat currents and the rectification are plotted as functions of the levels' mean energy . Panel (a), where the absolute values of the heat currents J + and J − are displayed, shows that both are bell-shaped with the maximum occurring near = 0, while they are strongly suppressed when | | increases beyond 5k B T. The bell-shape feature is a result of the compensation taking place in the open-circuit setup, whereby the charge current which would arise due to the temperature bias is counter-balanced by the appearance of a thermovoltage between the leads, effectively moving the weight of the curve towards the centre of the plot irrespective of the actual values of the temperature bias and of Δ . Consistently, the same behavior was found [56,57] in the thermal conductance of a multilevel interacting QD. Notice that the maximum occurs for a value of slightly away from zero. The height of the maximum, however, does depend on ΔT (linearly up to ΔT 1.25T) and the energy separation between the levels. Δ 2.8k B T is the value which yields the largest maximum for our choice of parameters. The height of the peak goes rapidly to zero by moving Δ away from this value, while the peak width remains virtually unaltered (in agreement with approximate analytical results for the thermal conductance in multilevel interacting QDs reported in references [56,57]). Such width is slightly increasing with ΔT, but only for ΔT > T, while it is virtually independent of tunneling constants. On the other hand, the rectification coefficient, plotted in panel (b) as a function of , shows a peculiar behavior: when > 5k B T the rectification goes to zero, while when < 5k B T the rectification tends to a finite value. The latter fact stems from the assumption that E C is the largest energy scale: indeed, as we will see in figure 11(b), a different behavior is found for finite E C . The behavior of heat currents and rectification with the bias ΔT is essentially the same as the one found in section 3.1.1.
Using similar arguments as for the closed-circuit setup, we could find that one can reach ideal rectification under the following conditions:

Cotunneling contributions
In the previous section we studied the rectification of a QD in the sequential tunneling regime and under the assumption that the charging energy E C is much larger than any other energy in the system. Such a condition allowed us to neglect the double occupation state of the QD and to find analytic expressions for the occupation probabilities (P 0 , P ↑ and P ↓ ) and currents (see equations (A5), (A8) and (A10)). In this section we include contributions from cotunnelling processes in the calculation of the current, allowing for a finite charging energy (i.e. for the double occupation of the QD). The latter account for coherent, second-order processes in the coupling Hamiltonian, that transfer an electron from one lead to the other via a virtual state either changing (inelastic) or not changing (elastic) the state of the QD. The cotunneling transition rates (for charge and energy) are calculated taking into account that the QD can be initially empty [(0, 0)], fully occupied [(1, 1)], or occupied by one electron [either (0, 1) or (1, 0)], where in the notation (i, j), i = 0, 1 refers to the level ↑ and j = 0, 1 refers to the level ↓ . Such cotunneling transition rates are calculated in details in appendix B. We stress that the inelastic cotunneling processes modify the state of the QD, thus modifying the MEs and their stationary solutions. Such modified MEs are reported in appendix B3, see equation (B32). Notice that the first square brackets on the right-hand-side of equations (B32) account for the sequential tunneling contribution only, where P 2 represents the probability for the QD to be doubly occupied.
The total currents are obtained by summing the cotunneling contributions to the sequential contribution. Since charge and energy currents are conserved, in the following we will focus only on the currents flowing out of the left lead and express them as and respectively. For consistency, in this section we account for the double occupancy of the QD even for the sequential currents. Therefore, I c seq (ΔT) and I seq (ΔT) are the currents calculated in the weak coupling regime, which, unlike equations (A8) and (A10), also accounts for the probability P 2 of finding two electrons in the QD and the related sequential processes. The expressions for the currents I c seq and I seq are reported in equations (B39) and (B40). Also the expressions for the cotunneling currents I c cot and I cot are collected in appendix B4.
Before discussing the results on the specific situations and setups, in the following we show that the heat current is symmetric with respect to = −E C /2. This can be understood by considering the symmetry properties of the Hamiltonian H QD . Indeed, equation (2) can be cast in the form where we have defined the operatorĥ σ = 1 −n σ . This proves that the Hamiltonian does not change by substituting σ with − σ − E C and replacing the operatorn σ with the operatorĥ σ . This means that the Hamiltonian H QD is particle-hole symmetric around −E C /2, implying that as long as E C is finite and provided that the average of the chemical potentials is zero (μ L + μ R = 0). Note that μ L = μ R = 0 in the closed-circuit setup and μ L = −μ R in the open-circuit setup.

Degenerate level
Let us first consider the case of a QD with a doubly-degenerate level, namely Δ = 0. A few observations are in order. The cotunnelling processes (either elastic or inelastic) in which both the initial and the final states have one electron in the same lead do not transfer energy, since Δ = 0, and can be ignored. Furthermore, the inelastic cotunneling processes that change the QD state from empty to doubly occupied occur rarely. Indeed, for such processes to happen, the proper initial conditions must be fulfilled, namely the leads must provide available electrons at high energy (order of E C ) and the QD must be empty (see the first line of equation (B46)). However, when the electrons in the leads have large enough energy, the QD is rarely empty because of the occurrence of sequential tunneling processes, while when the QD is empty, the electrons in the leads do not have enough energy to overcome E C . The same happens for the processes that empty the initially doubly occupied QD. Instead, for inelastic cotunneling processes that change the QD state from (1, 0) to (0, 1), and vice versa, the energy of the electrons involved in the process is enough to overcome the charging energy E C , making the process more likely to happen, as shown in figure 19. Therefore, the main contribution to the heat current comes from either elastic cotunneling processes or inelastic cotunneling processes that change the QD state (1, 0) ↔ (0, 1). However, at the end of this subsection we will find one scenario where cotunneling processes which move two electrons from/to the QD are responsible for the finiteness of the rectification.  correspond to the minima of the heat current in the sequential regime (dashed green line), while decrease the heat current for the values of which corresponds to the peaks.
Let us start discussing the sequential tunneling regime. The heat current is nearly zero at = 0 and at = −E C = −20k B T because the energy carried by the electrons which tunnel through the QD in these two cases is zero, as one can understand from equation (B40). Indeed, when = 0, the QD has vanishing probability to be doubly occupied [there is not enough thermal energy for the QD to be in the state (1, 1)], i.e. P 2 ≈ 0, while the functions F L↑ and F L↓ are also vanishing, since E C k B T. The remaining terms in equation (B40), however, account for tunneling of electrons which carry no heat because = 0. Similar arguments apply to the case = −E C . In this situation the QD is very likely occupied (namely The remaining terms account for tunneling of electrons which do not carry heat since + E C = 0. Cotunneling processes, however, allow electrons with energy different from to tunnel through the QD (through the virtual states), thus allowing a finite heat current to flow at = 0 and = −E C = −20k B T and giving rise to a reduction of the heat transfer. Far from resonance, when the sequential forward heat current J + seq decreases exponentially, the cotunneling processes become the dominant transport processes, increasing the heat current. Figure 7(a) also confirms that J + is symmetric with respect to The rectification coefficient is plotted in figure 7(b) as a function of , with a solid red curve when cotunneling contributions are included, and with a dashed green curve when sequential tunneling processes only are accounted for. First notice that the latter curve basically coincides with the curve in figure 3(c) in the range of energies considered. We note that cotunneling increases the rectification in the ranges of values of where it lowers the magnitude of the currents and decreases the rectification where it increases them. This behaviour can be understood as follows. Let us define the cotunneling corrections to the forward and backward heat currents as ΔJ ± = J ± − J ± seq , respectively. In turns out that the absolute values of ΔJ + and ΔJ − are nearly equal, but their sign is opposite. This happens because the main contribution to the heat currents, as notice above, comes from cotunneling processes that are elastic and from inelastic processes occurring when the QD is occupied by one electron. The cotunneling rates associated to such processes contain, under the integration symbol, the difference between the Fermi distributions of the leads (see for example equations (B8), (B26) and (B31)) and, therefore, change sign under the inversion of the temperature bias. Now, since ΔJ + ≈ −ΔJ − , we can express the effect of the cotunneling on the currents as J ± = J ± seq ± ΔJ + . This implies that the absolute value of both currents are either increased or decreased, depending on whether ΔJ + is positive or negative, respectively. Now, we can write the rectification with the cotunneling contributions as Therefore the rectification R coincides with the sequential tunneling regime's rectification R seq when the cotunneling correction is zero, namely ΔJ + = 0, is greater than R seq when ΔJ + < 0, and is smaller than R seq when ΔJ + > 0. In figure 8, the maxima over the QD level's energy of the rectifications with the cotunneling contributions R (solid red line) and without the cotunneling contributions R seq (dashed green line) are plotted as functions of the bias temperature ΔT. The dotted black line is the upper-bound of the rectification found in the sequential tunneling regime, see equation (11), which, in the case of Γ L = 2Γ R , is equal to 1/15 ≈ 0.067. We note that the cotunneling contributions increase the maximal rectification when the bias temperature is smaller than about 1.75T. The fact that for larger values of ΔT cotunneling contributions decrease R is a consequence of the fact that by increasing ΔT, the local maxima of both R and R seq , see figure 7(b), move towards = 0 and = −E C . For such values of , however, the cotunneling corrections ΔJ + to the heat currents are positive, see figure 7(a), thus producing, according to equation (19), a decrease of R with respect to R seq .
Open-circuit setup. When E C is finite, the heat current (in the sequential tunneling regime) can be finite in the open-circuit setup since charge current and energy currents are not proportional to each other (see equations (B39) and (B40)), contrary to what we found in section 3.1.1. This is due to the fact that, when E C is finite, the processes which involve the charging energy E C (the ones proportional to F ± Lσ in equation (B40)) transfer a different amount of energy, namely (E C + σ ), with respect to those which do not involve E C , while transferring the same charge. It is possible, however, to prove that the rectification still vanishes. Indeed, by substituting the solution of the master equations (B32), accounting for sequential tunneling processes only, in the expression of the charge current (equation (B39)), one finds that the condition which nullifies such current is f L F − R = f R F − L , independently of the tunneling constants Γ L and Γ R . By plugging in this condition into the expression of the energy current one finds that J + = −J − , i.e. there is no rectification. Cotunneling processes, however, can generate rectification.
In figure 9, forward heat currents (panel (a)) and rectification coefficient (panel (b)) are plotted as functions of the energy of the levels. In panel (a), the forward heat current (J + ), obtained including the cotunneling contributions, is plotted together with the heat current J + seq (= −J − seq ) relative to the sequential processes only. We first notice that J + seq is small but finite, as expected, although only around = −E C /2 = −10k B T, where the sequential processes involving the empty and doubly occupied states coexist. Indeed, around = −E C /2, P 0 and P 2 turn out to be both small but finite. Notice that for values of for which only one of the two (P 0 or P 2 ) is non-zero one finds that equations (B39) and (B40) are proportional to each other, implying that the heat current is zero. This is the case for < −E C = −20k B T, where one finds that P 0 ≈ 0 (the QD is at least singly occupied) and f − Lσ ≈ 0, or for > 0, where P 2 ≈ 0 (the QD cannot be doubly occupied) and F Lσ ≈ 0.
Remarkably, figure 9(a) shows that the cotunneling processes contribute substantially to the overall heat current, which then grossly deviates from the sequential result. On the one hand, this is due to the fact that cotunneling contributions to the heat current are essentially unrelated to the (overall) charge current, which is zero [58]. On the other hand, the cotunneling contributions to the charge current modify the thermovoltage (with respect to the sequential situation) that establishes between the leads in the open-circuit condition. In turn, such a thermovoltage influences the heat current (through the distribution functions entering the expressions of the sequential and cotunneling contributions). This leads to an additional indirect modification of the overall heat current, with respect to the sequential result. According to figure 9(a), J + reaches its maximum at = 0 and at = −E C = −20k B T, while being symmetric with respect to = −E C /2.
In figure 9(b), the rectification R which includes the cotunneling contributions is plotted as a function of the energy of the levels . R presents two maxima at | | ≈ 2.5k B T and at | − E C | ≈ 2.5k B T. However, the rectification is very small, at least one order of magnitude smaller than in the closed-circuit setup. The reason is that the main cotunneling contributions to the heat current, mentioned in the beginning of section 3.2.1, for degenerate levels change sign under the inversion of the temperature bias. More precisely, this is the case for the differences J u Lσ→Rσ − J u Rσ→Lσ in the second and third term of equation (B46), which are proportional to the energy integral of the Fermi functions, calculated at equal energy, of the two leads (see equations (B26) and (B31)). Moreover, it is important to mention that, for degenerate levels, P ↑ and P ↓ are virtually independent of the sign of the temperature bias. As already mentioned above, also the quantities J u ij,σ in equation (B44) change sign under the inversion of the temperature bias. In conclusion, rectification is generated by the rare cotunneling events, represented by the first and fourth term in equation (B46), that move two electrons from the lead to the QD or vice versa.

Non-degenerate levels
Let us now consider the case of a QD with two non-degenerate levels, namely Δ = 0. In this situation both elastic and inelastic co-tunnelling processes contribute significantly to the heat current.
Closed-circuit setup. In figure 10, the heat currents and rectification coefficient are plotted as functions of the average QD levels' energy . In panel (a), we plot the forward heat current which includes cotunneling contributions J + as a solid red curve and in the presence of sequential tunneling processes only J + seq as a dashed green curve. The latter curve closely resembles the one plotted in figure 4(a) for > −10k B T, meaning that double occupancy, at least for the value of E C considered, does not modify the results substantially. The behavior of the heat current is similar to the degenerate case, see figure 7, with a global minimum at = −E C /2 = −10k B T, and two symmetric local minima at ≈ 0 and ≈ −E C = −20k B T, which however do not touch zero, as in the degenerate case. We emphasize that the cotunneling contributions ΔJ + and −ΔJ − , unlike in the degenerate case, do not coincide (i.e. ΔJ + = −ΔJ − ). The reason for this is that the inelastic cotunneling processes that occur when the QD is occupied by one electron transfer a finite amount of heat and are not antisymmetric under the exchange of the leads' temperatures.
In figure 10(b) we plot the rectifications with (solid red curve, R) and without (dashed green curve, R seq ) cotunneling contributions as functions of the average QD levels' energy . The curve R seq , for > −10k B T, closely resembles the rectification reported in figure 4(b), calculated in the limit of infinite E C . Remarkably, the cotunneling contributions increase the rectification with respect to both R seq and the rectification  Open-circuit setup. In figure 11 we plot the heat currents and rectification as functions of the average QD levels' energy . In panel (a), we plot the forward heat current J + (solid red line) along with the one accounting only for sequential tunneling processes J + seq (dashed green line). The latter presents a local maximum at = −E C /2 = −10k B T of similar shape and height as in the degenerate case (see figure 9(a)). Notice that such a maximum does not appear in figure 6(a), where double occupancy of the QD was not allowed. J + seq peaks also at = 0, resembling the curve in figure 6(a), and = −E C = −20k B T: here heat transport is made possible by the energy difference Δ between the levels (i.e. charge and heat currents are not proportional to each other). Interestingly, the red curve is always above the green curve, meaning that cotunneling contributions increase the heat current for all values of . In addition, the main peaks are widened, while at = −E C /2 = −10k B T we have now a minimum. Also in this case the heat currents are symmetric around the axis = −E C /2 = −10k B T, as noticed in section 3.2.1.
In panel (b) of figure 11 we plot the rectification coefficients when cotunneling contributions are included (R, solid red curve) and when only sequential tunneling processes are allowed (R seq , dashed green curve). Notice that R seq virtually coincides with the rectification plotted in figure 6(b) (here plotted as a thin black curve) only for > −3k B T, whereas the two curves largely depart for other values of . This means that the rectification is much more sensitive to the finiteness of E C than the current, making clear that in the range −17k B T < < −3k B T the limit of infinite E C does not apply (electrons have the energy to overcome the Coulomb repulsion, represented by the value of E C , and the processes that involve E C can occur). As a result, R seq drops rapidly by decreasing below −3k B T, presenting a minimum at = −E C /2 = −10k B T, thus leaving a maximum at ≈ −4k B T.
As shown by the solid red curve, representing R, the cotunneling contributions affects very much the rectification by lowering it around the maxima of R seq and increasing it around = −E C /2 = −10k B T. In particular, R reaches its maximum at = −E C /2 = −10k B T, presents two symmetric local maxima at ≈ −2k B T and at ≈ −18k B T, and has two symmetric local minima at ≈ −4k B T and ≈ −16k B T.

Four-terminal device
The set up is shown in figures 1(c) and (d) and consists of a pair of Coulomb-coupled QDs, each attached to two leads. On the left-hand-side, QD ↑ is attached to L1, at temperature T + ΔT/2, and to L2, at temperature T − ΔT/2. On the right-hand-side, QD ↓ is attached to R1 and R2, both at temperature T. All reservoirs are kept at the same chemical potential, which is set to zero without loss of generality. There is natural flow of heat through QD ↑ from the hot to the cold reservoir, depending on the sign of ΔT. QD ↑ along with the two reservoirs attached to it constitute the drive circuit, while QD ↓ and reservoirs R1 and R2 constitute the drag circuit [59]. The drag circuit is coupled to the drive circuit via the Coulomb interaction: there is no particle exchange between the two circuits. The exchange of energy between the two circuits leads to a finite heat flow in the drag circuit [59]. Non-local heat rectification, i.e. rectification in the drag currents, occurs when the absolute value of the heat flowing between the drive and drag circuit depends on the sign of the temperature bias ΔT applied to the drive circuit. We must stress here that the currents J R1 and J R2 need not be equal, since there is an energy flow between drive and drag circuits. The sketch in figure 1(c) represents the forward bias configuration, with ΔT > 0 and drag currents indicated by J + R1 and J + R2 , while figure 1(d) represents the backward bias configuration, with ΔT < 0 and drag currents indicated by J − R1 and J − R2 . We fix the convention where heat currents are positive when they enter the QDs. The state of the system is described by the following set of occupancy (see section 2): (n ↑ , n ↓ ) = {00, 10, 01, 11}, where n σ represents the number of electrons in QD σ . Note that here, as in section 3.2, we allow for double occupation but we consider sequential tunneling only. As in appendix A and section 3.2, we describe the state of the system by the probabilities P 0 , P ↑ , P ↓ and P 2 , the latter referring to double occupancy. The MEs which allow to determine such probabilities are formally equal to equation (B32), reported in appendix B3. The heat currents (which coincides with the energy currents) relative to the drag circuit are given by where β = R1, R2. Before discussing the results, some general observations are in order. When, in the drive circuit, the coupling to the hot reservoir is stronger with respect to the coupling to the cold reservoir, and setting ↑ = ↓ , we notice that both currents in the drag circuit (J R1 and J R2 ) are negative (entering the leads), irrespective of all other parameters. In the opposite situation, where the coupling to the cold reservoir is stronger with respect to the coupling to the hot reservoir, the currents J R1 and J R2 are both positive (exiting the leads). This is not the case, however, when ↑ = ↓ , where the sign of J R1 and J R2 can be different and depend on all the parameters [60]. In particular, the sign of J R1 and J R2 does not depend on the direction of the temperature bias (forward or backward). However, when Γ R1 = Γ R2 , J R1 and J R2 are equal even when ↑ = ↓ and regardless of the values of Γ L1 and Γ L2 .
Remarkably, non-local rectification takes place only when the couplings in the drive circuit are asymmetric, i.e. when Γ L1 = Γ L2 . In what follows, for simplicity, we fix Γ R1 = Γ R2 = 0.05k B T, and we identify J ± ≡ J ± R1 = J ± R2 , with the rectification coefficient defined as in equation (1). In figure 12 we plot the heat current J (panel (a)) and the rectification coefficient (panel (b)) as a function of the energy level of the QD in the drive circuit ↑ in the presence of an asymmetry in the couplings in the drive circuit. For panel (a) the blue solid curve refers to the forward bias and the blue dashed-dotted curve refers to the backward bias, calculated at ΔT = ±0.3T, respectively. Both curves present a (positive) maximum around ↑ = 0 (the energy level is aligned with the chemical potential of the leads), in agreement with the results of reference [59]. Both currents vanish for large values of ↑ , since in this situation transport cannot occur even in the drive circuit, but for intermediate values of ↑ they show negative minima. Despite the relatively small difference between the coupling strengths in the drive circuit (Γ L1 = 0.05k B T and Γ L2 = 0.08k B T), the two curves depart significantly. This is quantified by the non-local rectification coefficient plotted in figure 12  of ↑ . J R1 and J R2 , however, are rather suppressed if compared with the heat currents relative to the setup with a single QD (figures 1(a) and (b)). Indeed, the heat current flowing in the drive circuit, in the case where drive and drag circuits are decoupled (E C = 0), turns out to be 3 orders of magnitude larger than in figure 12(a). A richer behavior occurs if one now assumes that the effective tunneling amplitudes t α (see equation (5)) are energy-dependent, thus depending on the charge state of the QDs. This situation was actually experimentally observed in reference [53], where the tunneling probabilities between the QD and the electron reservoirs in the drag circuit were found to depend on the charge state of the QD in the drive circuit. We can account for this situation by replacing the definition of the tunneling constants in the drag circuit (see equation (6)) with where the superscript (0) and (1) refer to the charge state (empty and occupied) of the QD in the drive circuit QD ↑ . The second equalities define the charge state-dependent coefficients κ  panel (b)). Notice that in this case J R1 and J R2 are different since Γ R1 = Γ R2 . Perfect heat rectification (R = 1) occurs also in this case for both currents J R1 and J R2 , and for the same values of ↑ (corresponding to points where the currents vanish). What is remarkable in figure 13 is that J R1 and J R2 have similar amplitude but opposite signs, both in the forward and backward temperature bias. This means that if heat is extracted from lead R1, a similar amount of heat is deposited into lead R2 (or the other way around). In particular, the extraction of heat from reservoir R1 can be used to lower its temperature, thus realizing an absorption refrigerator of the kind studied in reference [61], where cooling is driven by a non-local temperature difference, with no work provided to the system. We conclude the section by noting that we have checked that perfect rectification is obtained even when cotunneling processes are taken into account.

Conclusions
In this paper we have studied heat rectification in two different quantum dot (QD) systems: (i) a single QD in a two-terminal device and, (ii) a pair of coupled QDs in a four-terminal device. Heat rectification occurs when a QD is coupled asymmetrically to two terminals. Heat currents have been calculated using the master equation approach, up to the second order (cotunneling corrections) for a single QD. In case (i) we have considered a QD with either a doubly-degenerate level or two non-degenerate levels, each attached to two reservoirs. Furthermore, we have assumed that the device is either in the open-circuit setup, where the charge current vanishes, or in the closed-circuit setup, where the two reservoirs are kept at equal chemical potentials. In both cases energy current and heat current coincide. Within the sequential tunneling regime, we have first considered the case where the charging energy is very large such that the QD can only be occupied by a single electron. In this situation charge and heat current are proportional, so that no heat current can flow in the open-circuit setup in the case of a doubly-degenerate level. In the closed-circuit setup, however, the heat current is finite and, remarkably, we could derive an upper bound to the rectification coefficient R which only depends on the tunneling constants to the left and to the right as On the contrary, we have found that no bound exists in the case of two non-degenerate levels and we have identified the parameters' set which allow R to reach 1, its maximum value. Very interesting results are related to the effect of cotunneling contributions, including both elastic and inelastic processes, in the two-terminal configuration. In general, we have found that there exists ranges of values of the levels' position where rectification is enhanced by cotunneling. Moreover, we have found that • In the open-circuit setup of the degenerate level case, cotunneling processes, while permitting a finite heat flow, yield a finite (though small) rectification; • In the closed-circuit setup of the degenerate case, cotunneling corrections to the forward heat current are opposite to the corrections to the backward heat current, so that the magnitude of both currents are either increased or decreased by cotunneling. Rectification enhancement occurs when cotunneling lowers such magnitude (see also reference [33]); • In the open-circuit setup of the non-degenerate case, cotunneling always increases the currents with respect to the sequential regime. In the case (ii), we have considered the case of two Coulomb-coupled QDs, each connected to two terminals. This is a non-local configuration, where the temperature bias is applied to two terminals (the drive circuit), while the heat currents of interest are calculated in the other two terminals (the drag circuit). In this situation we have found, remarkably, that the rectification coefficient can reach the ideal value (i.e. R = 1), although with a rather small absolute value of the currents. Moreover, we have considered the experimentally relevant case [53] where the tunnel couplings between QDs and leads depend on energy. We have found that the heat currents in the drag circuit have similar amplitude but opposite signs, meaning that if heat is extracted from one reservoir, a similar amount of heat is deposited into the other.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix A. Master equation in the sequential tunneling regime
Here we assume that the charging energy E C is the largest energy scale, so that we can neglect all electronic configurations in which the total number of electrons in the QD exceeds one. We can describe the state of the QD by specifying the probability of finding the QD in the state with zero electrons P 0 , with one electron in a level with spin up (P ↑ ), and with one electron in a level with spin down (P ↓ ). The master equations can be written in matrix form as Here we have defined as the Fermi distribution function [f α (E)] of lead α = L, R evaluated at the QD levels' energies σ , while where The stationary master equation is obtained by equating the time derivative of P to zero and imposing the normalization of the probabilities P 0 + P ↑ + P ↓ = 1. We obtain the following stationary solutions Note that Λ changes non trivially under the exchange of the leads' temperatures. Such a behaviour results in a change in the QD state distribution which leads to rectification. We have also calculated the master equation accounting for up to two electrons in the QD. When E C is two order of magnitude larger than k B T we have proven that the probability of occupation of the states with two electrons is negligible.
The charge current entering the QD from the left lead can be written as (−e is the electronic charge) and becomes after substituting the solutions of the master equation, equation (A5). Similarly, the energy current takes the form and becomes Finally, from the definition of heat current

Appendix B. Cotunneling contributions for a single QD with two states
We consider the QD in the weak coupling condition, so that we can treat the tunnel Hamiltonian H T in equation (4) as a perturbation to the system. Therefore, we describe the system through the eigenstates of the free Hamiltonian H 0 = H QD + H L + H R and calculate the transition rates between two of such states using the generalized Fermi golden rule where γ if is the rate associated with the process that starts from the initial state |i with energy E i , and arrives in the final state |f with energy E f . Since the perturbation H T is time-independent, the delta function in equation (B1) imposes the energy conservation between the initial and final states E i = E f . Moreover, the amplitude A if contains the perturbation term H T which gives rise to a natural expansion in the powers of H T . The first order describes the sequential tunneling processes, while the second order describes the cotunneling processes [63,64]. In this appendix we calculate the transition amplitudes of the latter processes.
In a cotunneling process, the system evolves from an initial state to a final one passing through a virtual state. Since there can be more than one virtual state, the cotunneling processes can exhibit quantum interference. The transition amplitude A if that enters the generalized Fermi golden rule can be written as where the sum is made over all the virtual states of the system, and the parameter η goes to zero and is needed to eliminate the divergences in the calculation of the rates [65], see appendix C. Such divergences are due to the sequential tunneling processes. Indeed, the system can evolve from the initial state to the final state also through two consecutive sequential tunneling events. For example, an electron of the left lead can tunnel sequentially into the QD and, from there, it can tunnel sequentially into the right lead. Such a process transfers an electron from the left to the right one, but it is made of two sequential tunneling events. When integrating over all possible initial and final states, both the transition rates that are given by the second-order Fermi golden rule and the processes made of two consecutive sequential tunneling events contribute to the integral. However, since we have already accounted for the sequential tunneling, we have to remove the contributions of the sequential tunneling events from the transition rates and the currents of the cotunneling removing the divergences of the integrals. For every pair of initial and final states, we have to calculate the transition rate of each given process and the currents associated with it by multiplying the Fermi golden rule rate by the Fermi distributions relevant to the tunneling process. The Fermi distributions are necessary to describe the probability of having the starting electrons in the initial state and the final electron levels empty so that they can be occupied by the incoming electrons. Then, we obtain the total rates by summing over all possible initial and final states.
We separate the cotunneling processes in two kinds: the elastic processes, in which the energy of the QD does not change between the initial and final state, and the inelastic processes, which modify the QD energy. Thus, the inelastic processes modify the state of the QD and enter the master equation. Whereas, the elastic processes do not. However, both elastic and inelastic processes contribute to the transport of charge and heat.

B1. Cotunnelling rates: elastic processes
In this section we derive the cotunneling rates and the currents for the various elastic processes. The elastic processes do not change the energy of the QD, therefore the initial and final states of the QD must be the same. Since the QD can be occupied by either zero, one, or two electrons, we separate the cotunneling contributions according to the QD state. For each QD state, we find the possible initial and final states and the corresponding transition rate. Of course, we do not consider the processes in which the initial and final states are in the same lead because such processes do not contribute to the transport of heat nor charge. In general, the cotunneling rate for an electron to go from lead α to lead β, while the QD is initially in the state 'in', can be calculated as where D α is the density of states of lead α and the Fermi distribution f α describes the probability of finding an occupied electronic state in the lead α, while f − β = 1 − f β is the probability of finding an unoccupied state in the lead β. Of course, the cotunneling rate for the opposite process, i.e. for an electron to go from lead β to lead α, is obtained by exchanging the leads indices in equation (B4).
On the other hand, to calculate a net current we have to account for both the L → R and the R → L processes. The net single-process charge current, from left to right, can thus be written as while the net single-process energy current, from left to right, can be written as where e is the electron charge and 'in' refers to the state of the QD (with spin σ).
In the following we list the expressions for the transition rates and the currents depending on the initial state of the QD. We will use the superscript ij (with i, j = 0, 1) to indicate that the QD is initially in the configuration (in) =(i, j), i.e. there are i electrons in the level ↑ and j electrons in the level ↓. Namely, when  the QD is initially empty and in its intermediate state the level σ is occupied, see for example figure 14, we find and where K = −e for the charge current (superscript c) and K = E for the energy current (superscript u). When the QD is initially fully occupied and in the intermediate state the level σ is empty, we find and When the QD is initially occupied by one electron with spin up and in its intermediate state the level σ is empty we find and Finally, when the QD is initially occupied by one electron with spin down and in its intermediate state the level σ is empty, see for example figure 15, we find and After removing the divergent part of the integrals of equations (B7), (B9), (B11) and (B13), see appendix C for the details, the transition rates of the cotunneling processes can be negative. Despite this, the total transition rate, which also accounts for two consecutive sequential tunneling events, is always positive as it must be.

B2. Cotunnelling rates: inelastic processes
In this section we derive the cotunneling rates and the currents for the various inelastic processes. In an inelastic process the state of the QD gets modified, so that the energy of the QD changes. This can take place in two different ways only, either by adding/removing two electrons to the QD, or by removing an electron from one level and adding one in the other level. In particular, the change of states are • (0, 0) → (1, 1), the QD is initially empty and, through the inelastic process, becomes fully occupied; • (1, 1) → (0, 0), the QD is initially fully occupied and, through the inelastic process, becomes empty; • (1, 0) ↔ (0, 1), in both initial and final state the QD has one electron inside, but the inelastic process changes the occupied level. Analogously to appendix B1, we can organize the different processes on the basis of the initial state of the QD.
When the QD is initially empty, the only possible final state is the one in which there are two electrons in the QD, since cotunneling comprises two tunneling processes. The two electrons can both come from the left lead (i), or from the right lead (ii), or one from the left and one from the right lead (iii). In case (i), there are two possible intermediate states, depending on which level is occupied first in the QD, as shown in figure 16. Imposing the energy conservation between the initial and the final states, we obtain the following cotunneling rate Since the electrons involved in the process carry both charge and energy out of the left lead, the current takes the form where K = −2e for the charge current (superscript c), since both electrons that tunnel into the QD come from the left lead, and K = ↑ + ↓ + E C for the energy current (superscript u), which is the energy removed from the left lead in the cotunneling process. Notice that K does not depend on the integration variable, therefore the charge and energy currents are proportional to each other. Case (ii) is analogous to case (i) and we can calculate the cotunneling rate by exchanging the leads' insides L → R in equation (B15) obtaining Note that there is no current associated to case (ii), since neither electron tunnels from the left lead.
In case (iii) there are two possible initial states, depending on whether the electron with spin up comes from the left or right lead. For each initial state there are two possible intermediate states, depending on whether the electron occupying the QD comes from the left or the right lead. Imposing the energy conservation between the initial and the final states, we obtain the following cotunneling rate where K = −e for the charge current (superscript c) and K = E for the energy current (superscript u).
Notice that in this case charge and energy currents are not proportional to each other since K = E cannot be taken out of the integration. Let us now consider the case where the QD is initially occupied by two electrons. The inelastic cotunneling processes, that move both electrons out of the QD, are the inverse processes with respect to the one discussed above (relative to the QD initially empty). Therefore, the initial and the final state of the processes which empty the QD are, respectively, the final and the initial state of the processes which fill the QD. Moreover, the intermediate states are the same. Therefore, we obtain the following expressions for the cotunneling rates which is analogous to the transition rate of equation (B18), and For the currents flowing out of the left lead we obtain where K = −2e for the charge current (superscript c) and K = ↑ + ↓ + E C for the energy current (superscript u), where K = −e for the charge current (superscript c) and K = E for the energy current (superscript u). Finally, let us now consider the case where initially a certain level of the QD is occupied by one electron. After the inelastic processes, in the final state the QD will still contain one electron, but in other level. When one of the electron is initially in the left lead, there are two possible initial states, one for each value of the spin. For each initial state there are two final states: (i) the electron in the QD tunnels into the right lead, so that both charge and energy are transferred in the cotunneling process; (ii) the electron in the QD tunnels into the left lead, in which case only energy is transferred (the QD levels have different energies). In case (i), see figure 17, there are two possible intermediate states (QD fully occupied and QD empty), depending on the order of the two tunneling processes. We obtain the following cotunneling rates These processes transfer both charge and energy, so that the current flowing out of the left lead is given by where K = −e for the charge current (superscript c) and K = E for the energy current (superscript u). Also in case (ii), see figure 18, there are two possible intermediate states (QD fully occupied and QD empty), depending on the order of the two tunneling processes. We obtain the following cotunneling rates (B27) Figure 17. Diagram of the inelastic cotunneling process L → R with incoming electron with spin ↑ and QD initially in (0, 1). These processes transfer only energy, with a current flowing out of the left lead given by where K = σ − σ . When one of the electron is initially in the right lead, analogous calculations lead to the following cotunneling rates and expressions of the current flowing out of the left lead (only for the former processes, since the latter processes do not produce any current in the left lead) where K = −e for the charge current (superscript c) and K = E for the energy current (superscript u).
Finally we notice that, in the present case, cotunneling inelastic processes occur via two intermediate states (contrary to elastic processes, which occur via a single intermediate state), thus giving rise to quantum interference effects.

B3. Master equations and cotunnelling rates
When inelastic co-tunnelling processes are included, the MEs need to be modified and can be written as follows d dt where is the sum of the inelastic cotunneling transition rates relative to the processes that empty the QD, and is the sum of the inelastic cotunneling rates of the processes that fill the QD. Similarly, T 01→10 and T 10→01 are defined as the sums of the inelastic cotunneling rates that exchange the level in the QD which is occupied.
In equation (B32), P 0 represents the probability for the QD to be unoccupied, P 2 represents the probability for the QD to be doubly occupied, while P ↑ (P ↓ ) is the probability for the lower ↑ (upper ↓ ) level of the QD to be occupied. Moreover, we have defined and where Θ + σ describes the total rate of the sequential processes that move an electron with spin σ from the leads to the QD that has already one electron inside, while Θ − σ represents the inverse process, in which the electron with spin σ tunnels from the doubly occupied QD to the leads. Notice that, in equation (B32), the first square brackets contain the sequential tunneling contribution, see equation (A1), whereas, the second square brackets refer to the cotunneling corrections. The equation relative to P ↓ is obtained by exchanging the labels ↑ and ↓ in the last line of equation (B32).
In the four-terminal device, the MEs which allow to determine the probabilities P 0 , P ↑ , P ↓ and P 2 in the sequential tunneling regime are formally equal to equation (B32), where we consider only the first square brackets in each line. In this case, however, the quantities Θ ± σ and Σ ± σ are defined as follows and and the tunneling constants Γ α are defined by equation (6), with α = L1, L2, R1, R2.

B4. Expressions for the currents
In the case of a QD with two levels one finds and where F Lσ = f L ( σ + E C ) and F − Lσ = 1 − f L ( σ + E C ). Notice that equations (B39) and (B40) reduce to equations (A7) and (A9), respectively, when E C diverges, thus causing F Lσ and P 2 to vanish.
In equation (17) while the inelastic components are and Such probabilities are calculated through a ME which also account for inelastic co-tunnelling processes (see appendix B3). The elastic co-tunnelling single-process currents (J It is worth noticing that the elastic cotunneling contributions to the heat current can give rise to rectification, despite the fact that the quantities J ij,u σ in equation (B44) depend (under an energy integration) on the difference between the Fermi functions of the two leads at the same energy. Indeed, the probabilities P 0 , P 2 , P ↑ and P ↓ actually depend on the sign of the temperature bias, in a more noticeable way for non-degenerate levels.
Furthermore, we can identify each line of equation (B46) with the contribution to the energy current of the corresponding change of QD state, namely The contributions J 00→11 and J 11→00 are suppressed because the respective single-process currents and QD occupation probabilities decrease exponentially in different energy regions. Indeed, for the process (0, 0) → (1, 1) [(1, 1) → (0, 0)], the probability P 0 (P 2 ) decreases as the electrons (holes) from the leads can occupy the QD levels, just through sequential tunneling, when 0 ( −E C ). On the other hand, the single-process currents decreases as the leads electrons (holes) do not have enough energy to overcome the charging energy E C when −E C /2 ( −E C /2). Instead, the cotunneling processes (1, 0) ↔ (0, 1) are not much suppressed because the combined energy of the involved electrons does not have to overcome the charging energy E C . In figure 19 these inelastic cotunneling contributions to the heat currents are compared: (1, 0) ↔ (0, 1) and (0, 0) ↔ (1, 1) contributions differ in order of magnitudes.

Appendix C. Renormalization of cotunneling integrals
To calculate the transition rates and the currents of the cotunneling processes, we have to integrate over all the initial and final states. We note that all these integrals diverge when η → 0, so we have to separate and remove the divergent part from the rest. As shown below, such divergent part is associated with the sequential tunneling [65]. A generic integral is in the following form: where A and B are some constant energies and g(E) is a smooth function of E. Expanding the squared modulus, we can separate the integral in three parts g(E)dE, where the integrals I 1 contain the squared moduli of the fractions, while the integral I 2 contains the real part of the product between the fractions. Let us see how to regularize the integral I 1 . Expanding the squared modulus and changing the integration variable, we obtain Then, we can sum and subtract the quantity g(A) at the numerator of the integral, obtaining The first integral can be calculated exactly and is proportional to 1/η, namely The first term is divergent in the limit η → 0, so we remove it. Whereas the second term converges at vanishing η. Then, the integral I 1 is reduced to The integrand diverges at E = 0, so, we need to eliminate this divergence to compute numerically the integrals. We note that, changing the sign of the integration variable, it yields Therefore, summing the integrals, we obtain in which the integrand does not diverge at E = 0 and can be computed numerically. We note that after we remove the divergence proportional to 1/η, the sign of the integral I 1 (A) is no longer guaranteed to be positive. Therefore, the transition rates of the cotunneling processes can have negative values. However, the total transition rate of the processes that move the system from the initial state to the final state, namely the sequential tunneling processes and the cotunneling processes, is always positive, as it is given by equation (C5). Now, let us regularize the integral I 2 . First, we expand the real part of the product of the fractions, obtaining In the limit η → 0, the term proportional to (E − A)(E − B) becomes the principal values of the integral, while the term proportional to η 2 vanishes. Therefore, the integral I 2 is reduced to where the integral is performed on the principal value. To compute the integral numerically, we can split the denominator in the following way: Then, we can split the integral and change variables to move the singularity in E = 0, namely Finally, we can change the integration variable from E to −E and obtain a similar integral with opposite sign. Summing these integrals together, we obtain which is computable numerically because the integrand has no divergences. We remove the divergent part because it is generated by the sequential processes. Indeed, in the term proportional to 1/η there is energy conservation between the initial, the intermediate, and the final state. Therefore, the system can arrive in the final state through two sequential tunneling processes. Moreover, the imaginary parameter iη is associated with the rate of leaving the intermediate state [66]. Indeed, imaginary energies describe metastable states and decay processes. For instance, suppose to initialize the system in a state |α , which has energy E − iη. Then, after a time t, the state has evolved into |α t = e −i(E−iη)t/ |α , and the probability of finding the system in the state |α becomes e −2ηt/ . In the same way, η is due to the processes that move the system out of the state |α . For example, let us consider a QD with two non-degenerate levels of energies ± Δ /2 and let us calculate the divergent part of the cotunneling processes that have as the intermediate state the state in which the QD has one electron of spin σ. Using the currents calculated in appendix B1, and considering the divergent part that arises from equation (C5), we can calculate the divergent part of the elastic cotunneling charge current that leaves the left lead, namely Then, we obtain the divergence of the inelastic cotunneling processes using the currents calculated in appendix B2, namely Summing together such currents, we get Now, we remember that the rate equation of the probability P σ in the sequential tunneling regime is and the left charge current transferred by the processes that enter such a rate equation is In the stationary condition, the charge current becomes Comparing equations (C17) and (C20), we can find that η satisfies which is half the rate of the sequential tunneling processes that move the QD out of the state with one electron of spin σ. This state is the intermediate state of the cotunneling processes considered.