Thermoelectrics in Coulomb-coupled quantum dots: Cotunneling and energy-dependent lead couplings
TThermoelectrics in Coulomb-coupled quantum dots:Cotunneling and energy-dependent lead couplings
Nicklas Walldorf, ∗ Antti-Pekka Jauho, and Kristen Kaasbjerg † Center for Nanostructured Graphene (CNG), Dept. of Micro- and Nanotechnology,Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark (Dated: September 20, 2018)We study thermoelectric effects in Coulomb-coupled quantum-dot (CCQD) systems beyondlowest-order tunneling processes and the often applied wide-band approximation. To this end,we present a master-equation (ME) approach based on a perturbative T -matrix calculation of thecharge and heat tunneling rates and transport currents. Applying the method to transport througha non-interacting single-level QD, we demonstrate excellent agreement with the Landauer-B¨uttikertheory when higher-order (cotunneling) processes are included in the ME. Next, we study the ef-fect of cotunneling and energy-dependent lead couplings on the heat currents in a system of twoCoulomb-coupled QDs. Overall, we find that cotunneling processes (i) dominate the heat currentsat low temperature and bias, and (ii) give rise to a pronounced reduction of the cooling powerachievable with the recently demonstrated Maxwell’s demon cooling mechanism. Furthermore, wedemonstrate that the cooling power can be boosted significantly by carefully engineering the energydependence of the lead couplings to filter out undesired transport processes. Our findings empha-size the importance of considering higher-order cotunneling processes as well as the advantage ofengineered energy-dependent lead couplings in the optimization of the thermoelectric performanceof Coulomb-coupled QD systems. I. INTRODUCTION
The experimental progress in control of single-electrontransport has spurred interest in nanosystems which uti-lize the associated heat currents for thermoelectric ap-plications. In particular, experiments with Coulomb-coupled quantum-dot (CCQD) systems have demon-strated a plethora of novel phenomena ranging fromCoulomb drag and electron pairing to extraordinarythermoelectric effects. This includes the realization ofan energy harvester which converts a thermal gradientin a CCQD system into an electric current, as well asan autonomous Maxwell’s demon capable of cooling acurrent-carrying QD system at the cost of heating a “de-mon” QD system. In addition to the above, theoretical studies have pre-dicted a wide range of novel thermoelectric effects inCCQD systems.
The mechanisms behind these ef-fects rely on the presence of a strong Coulomb interac-tion between electrons in the otherwise decoupled QDs(see Fig. 1 for the case of two Coulomb-coupled QDs).The strong interaction can be utilized to tailor the ther-moelectric properties of CCQD systems and opens theopportunity to test fundamental thermodynamic aspectsof heat transport in interacting nanoscale systems drivenout of equilibrium. While the operation principles of the above-mentionedeffects are governed by incoherent electron tunneling (se-quential tunneling) processes between the leads and theQDs, the importance of coherent higher-order tunnel-ing (cotunneling) processes for the nonlinear heat trans-port remains largely unexplored. Furthermore, when op-erated under strong non-equilibrium conditions wherelinear response theory breaks down, a theoretical treat-ment taking into account the full nonlinear properties is needed.
Only recently have these issues been dis-cussed in strongly interacting QD systems.
Another important factor for thermoelectric effects inCCQD systems is the coupling to the leads which isusually treated in the wide-band approximation assum-ing energy-independent couplings. However, energy-dependent couplings to the leads occur naturally in manyQD systems, and add an important degree of tun-ability to the system, and is as crucial for the thermoelec-tric properties as it is for Coulomb drag.
In this work, we present a master-equation approachfor the calculation of the nonlinear electronic charge andheat currents in interacting QD systems which takes intoaccount the above-mentioned factors. The charge andheat transfer rates produced by electron tunneling pro-cesses are obtained with a perturbative T -matrix ap-proach which allows us to treat sequential and cotun-neling processes on equal footing. We resolve the techni-cal challenges associated with the evaluation of the cotun-neling rates with an implementation of the often appliedregularization scheme which applies to the generalcase of energy-dependent lead couplings, applied biases,and temperature gradients in the system.The main findings and the organization of the paperare as follows. In Sec. II, we introduce the model sys-tem of CCQDs. In Sec. III, we present the methodol-ogy, and benchmark the approach in Sec. IV by com-paring it to the Landauer-B¨uttiker formalism for trans-port through a non-interacting single-level QD. In Sec.V, we study nonlinear thermoelectric phenomena in CC-QDs. We investigate the energy exchange mediated bythe inter-dot Coulomb interaction which among otherthermoelectric effects leads to the demon-induced coolingmechanism. Our findings shed light on the limitationsimposed by cotunneling processes on the performance of a r X i v : . [ c ond - m a t . m e s - h a ll ] J un this mechanism. Furthermore, we demonstrate a stronglyenhanced performance of the demon-induced cooling ef-fect by tuning the energy-dependence of the lead cou-plings. In such performance optimization, as we show,cotunneling processes are essential for a quantitative de-scription of the thermoelectric properties. Finally, Sec.VI presents our conclusions, and App. A gives technicaldetails on the cotunneling rates and the regularizationprocedure. II. COULOMB-COUPLED QD SYSTEMS
We consider CCQD systems like the one illustrated inFig. 1, which can be described by the Hamiltonianˆ H = ˆ H dots + ˆ H leads + ˆ H T , (1)and consists of a system of CCQDs with Hamiltonianˆ H dots which is coupled to external leads with Hamilto-nian ˆ H leads by tunnel couplings described by ˆ H T . Wedenote ˆ H = ˆ H dots + ˆ H leads .For the QD system, we consider a minimal spinlessmodel of inter-dot Coulomb-coupled single-level QDs de-scribed by the Hamiltonianˆ H dots = (cid:88) δ (cid:15) δ ˆ c † δ ˆ c δ + (cid:88) (cid:104) δ,δ (cid:48) (cid:105) U δδ (cid:48) ˆ n δ ˆ n δ (cid:48) , (2)where ˆ c † δ (ˆ c δ ) creates (annihilates) an electron in QD δ with energy controlled by gate voltages (cid:15) δ = − eV δ , where V δ is the gate potential on dot δ , ˆ n δ = ˆ c † δ ˆ c δ is the occu-pation number operator, U δδ (cid:48) is the inter-dot Coulombinteraction, and the summation in the second term isover all QD pairs (specific systems are studied in Secs.IV–V). U µ A T B µ B T C , µ C System 1System 2 T A A B
QD1 C QD2
FIG. 1. Illustration of the CCQD system studied in Sec. Vconsisting of two Coulomb-coupled QDs δ ∈ { , } withinter-dot Coulomb interaction U , tunnel-coupled in a three-terminal configuration to leads (cid:96) ∈ { A, B, C } (no tunnelingallowed between the QDs) with temperatures T (cid:96) and electro-chemical potentials µ (cid:96) . The leads are described by non-interacting electronreservoirs, ˆ H leads = (cid:80) (cid:96)k (cid:15) (cid:96)k ˆ c † (cid:96)k ˆ c (cid:96)k , where ˆ c † (cid:96)k (ˆ c (cid:96)k ) creates(annihilates) an electron with momentum k and energy (cid:15) (cid:96)k in lead (cid:96) , which is assumed to be in local equilib-rium with temperature T (cid:96) and electrochemical potential µ (cid:96) = µ − eV (cid:96) , where µ is the equilibrium chemical po-tential and V (cid:96) is the voltage applied to lead (cid:96) . The tun-neling Hamiltonian which couples the QD system to theleads is ˆ H T = (cid:80) (cid:96)kδ ( t (cid:96)kδ ˆ c † δ ˆ c (cid:96)k + h.c. ) , where t (cid:96)kδ is thetunneling amplitude. We define lead coupling strengthsas γ (cid:96) ( (cid:15) ) ≡ πd (cid:96) ( (cid:15) ) | t (cid:96) ( (cid:15) ) | , where d (cid:96) ( (cid:15) ) is the lead den-sity of states. γ (cid:96) ( (cid:15) ) is allowed to be energy dependent incontrast to the often applied wide-band approximation. III. MASTER EQUATION AND TRANSPORTCURRENTS
We describe the dynamics and transport in the CCQDsystem with a Pauli ME where the transitions betweenthe QD states are governed by electron tunneling toand from the leads. The tunneling-induced transitionrates are calculated based on a perturbative T -matrixapproach where the tunneling Hamiltonian is treated asa perturbation to the decoupled QD system and leads.This allows a systematic expansion in the tunnel cou-plings and the inclusion of high-order processes. How-ever, quantum effects such as tunneling-induced levelbroadening and level shifts are not captured by thisperturbative approach, which is only valid in the weakcoupling regime γ < k B T, U .In the absence of tunnel coupling, the states of the de-coupled QD system and leads are described by productstates of the QD system occupation states | m (cid:105) with en-ergy E dots ,m = (cid:104) m | ˆ H dots | m (cid:105) and the leads | i (cid:105) with energy E leads ,i = (cid:104) i | ˆ H leads | i (cid:105) . The non-equilibrium occupationsof the QD states are described by probabilities p m (thediagonal components of the reduced density operator ofthe CCQD system) which are determined by the ME˙ p m = (cid:88) n (cid:54) = m (Γ nm p n − Γ mn p m ) , (cid:88) m p m = 1 , (3)where Γ mn denotes the tunneling-induced transition ratefrom QD state | m (cid:105) to | n (cid:105) . The ME is solved for thesteady-state probabilities, ˙ p m = 0, in the following. TheQD states are given explicitly in Sec. IV and Sec. V forthe considered systems. A. Transition rates
The rates for transitions between the QD states areobtained from the generalized Fermi’s golden rule ˜Γ mn = 2 π (cid:126) (cid:88) ij |(cid:104) j |(cid:104) n | T | m (cid:105)| i (cid:105)| ρ i δ (∆ mn + E leads ,j − E leads ,i ) , (4)where ∆ mn ≡ E dots ,n − E dots ,m , ρ i is the thermal proba-bility of finding the leads in the initial state, the sum isover initial and final states of the leads, and the T matrixobeys ˆ T = ˆ H T + ˆ H T E initial − ˆ H + iη ˆ T , (5)with E initial = E dots ,m + E leads ,i , and η is a positive in-finitesimal.The lowest-order contribution to the tunneling ratesdescribes single-electron tunneling, or sequential tunnel-ing , processes between the QD system and the leads:Γ −→ (cid:96) mn = (cid:126) − γ (cid:96) (∆ mn ) f (cid:96) (∆ mn ) , (6)Γ ←− (cid:96) mn = (cid:126) − γ (cid:96) (∆ nm ) ¯ f (cid:96) (∆ nm ) , (7)where Eq. (6) (Eq. (7)) is the sequential rateof tunneling out of, → , (into, ← ) lead (cid:96) , therebychanging the state of the QD system from m to n , f (cid:96) ( (cid:15) ) = [exp ( β (cid:96) ( (cid:15) − µ (cid:96) )) + 1] − is the Fermi-Dirac distri-bution in lead (cid:96) , ¯ f (cid:96) ( (cid:15) ) = 1 − f (cid:96) ( (cid:15) ), and β (cid:96) = 1 / ( k B T (cid:96) ).The leads are assumed to equilibrate to the Fermi-Diracdistribution in between the tunneling events.The next-to-leading order terms in the T matrix de-scribe cotunneling processes. In conventional local elas-tic and inelastic cotunneling processes, a net electron istransferred between two leads attached to the same QD(e.g., System 1 in Fig. 1). Here we also consider (i) non-local cotunneling processes in which a net electron istransferred between leads attached to different QDs, aswell as (ii) pair-cotunneling processes where two electronstunnel into/out of the CCQD system in one coherent pro-cess. For the thermoelectric effects in focus here, the pro-cess of nonlocal cotunneling is important. The (unregu-larized) rate for nonlocal cotunneling which net transfersan electron out of lead (cid:96) and into lead (cid:96) (cid:48) is given by˜Γ −→ (cid:96) ←− (cid:96) (cid:48) mn = (cid:90) d(cid:15) π (cid:126) γ (cid:96) ( (cid:15) ) γ (cid:96) (cid:48) ( (cid:15) − ∆ mn ) f (cid:96) ( (cid:15) ) ¯ f (cid:96) (cid:48) ( (cid:15) − ∆ mn ) × (cid:12)(cid:12)(cid:12)(cid:12) vm + (cid:15) + iη + 1∆ v (cid:48) n − (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (8)where v ( v (cid:48) ) refers to the virtually occupied intermediatestate in the process where an electron initially tunnelsfrom lead (cid:96) and into the QD system (from the QD systemand into lead (cid:96) (cid:48) ). We refer to App. A for the expressionsfor the remaining cotunneling processes relevant for thisstudy.A well-known artifact of the cotunneling rates obtainedwith the T -matrix approach is that they formally divergein the limit η →
0. To deal with this divergence differ-ent regularization schemes have been proposed.
Deep inside the Coulomb blockade, the discripancy be-tween the different regularization schemes vanishes. Inthis work, we apply the by now standard regularizationscheme in Ref. 29, but for future work, a detailed compar-ison of the charge and heat currents obtained from dif-ferent regularization schemes could be useful. We denote the regularized rates which enter into Eq. (3) without atilde. To be explicit, we consider the processes Γ mn ≡ (cid:80) (cid:96) (Γ (cid:96) ← mn +Γ (cid:96) → mn ) , Γ (cid:96) ← mn ≡ Γ ←− (cid:96) mn + (cid:80) (cid:96) (cid:48) (Γ ←− (cid:96) −→ (cid:96) (cid:48) mn +Γ ←− (cid:96) ←− (cid:96) (cid:48) mn ) , Γ (cid:96) → mn ≡ Γ −→ (cid:96) mn + (cid:80) (cid:96) (cid:48) (Γ −→ (cid:96) ←− (cid:96) (cid:48) mn + Γ −→ (cid:96) −→ (cid:96) (cid:48) mn ) . A numerical procedure for theregularization is outlined in App. A.
B. Charge and heat currents
The steady-state transport currents can be obtainedfrom the occupation probabilities. The electric currentgoing into lead (cid:96) is I (cid:96) ≡ − e (cid:42)(cid:88) k d ˆ n (cid:96)k dt (cid:43) = − e (cid:88) mn p m (cid:0) Γ (cid:96) ← mn − Γ (cid:96) → mn (cid:1) , (9)where ˆ n (cid:96)k = ˆ c † (cid:96)k ˆ c (cid:96)k , p m is calculated from Eq. (3), andthe rightmost form expresses the electric current in termsof the total rate of electrons tunneling into lead (cid:96) , minusthe total rate of electrons tunneling out of lead (cid:96) . The heat current going into lead (cid:96) is J (cid:96) ≡ (cid:42)(cid:88) k ( (cid:15) (cid:96)k − µ (cid:96) ) d ˆ n (cid:96)k dt (cid:43) = (cid:88) mn p m (cid:0) W (cid:96) ← mn − W (cid:96) → mn (cid:1) , (10)where the rightmost form expresses the heat current interms of heat rates W (using a similar notation as for thetunneling rates).The sequential-tunneling heat rate in lead (cid:96) is calcu-lated as the tunneling rate multiplied by the energy ofthe tunneling electron relative to the chemical potentialin the lead, W −→ (cid:96) (cid:96),mn = (∆ mn − µ (cid:96) )Γ −→ (cid:96) mn ,W ←− (cid:96) (cid:96),mn = (∆ nm − µ (cid:96) )Γ ←− (cid:96) mn , (11)where the indices follow the notation of the tunnelingrates, however, the additional first subscript (cid:96) refers tothe lead in which the heat rate is calculated.Analogously, the cotunneling heat rates into/out of theleads are calculated a posteriori by multiplying the inte-grand in the cotunneling rate by the energy of the tun-neling electron relative to the chemical potential of thelead. For example, for the nonlocal cotunneling processbetween lead (cid:96) and (cid:96) (cid:48) , the heat rate in lead (cid:96) reads˜ W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96),mn = (cid:90) d(cid:15) π (cid:126) γ (cid:96) ( (cid:15) ) γ (cid:96) (cid:48) ( (cid:15) − ∆ mn ) f (cid:96) ( (cid:15) ) ¯ f (cid:96) (cid:48) ( (cid:15) − ∆ mn ) × ( (cid:15) − µ (cid:96) ) (cid:12)(cid:12)(cid:12)(cid:12) vm + (cid:15) + iη + 1∆ v (cid:48) n − (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (12)with the heat rate in lead (cid:96) (cid:48) , ˜ W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96) (cid:48) ,mn , given as above butwith ( (cid:15) − µ (cid:96) ) replaced by ( (cid:15) − ∆ mn − µ (cid:96) (cid:48) ). The remainingcotunneling heat rates follow similarly.Whereas the calculation of charge currents involves theelectron-tunneling rates which enter the ME (3), andtherefore does not require any additional steps once theME has been set up and solved, the heat currents mustbe calculated via the heat tunneling rates in a post-processing step, similar to more rigorous density-matrixtreatments. IV. COMPARISON TO THELANDAUER-B ¨UTTIKER FORMALISM
In this section, we benchmark the approach by com-paring the charge and heat currents in a spinless non-interacting single-level QD system with those obtainedfrom the Landauer-B¨uttiker (LB) formalism (see Ref. 42for a comparison of the electric current in the case ofequal temperatures in the leads). For non-interactingsystems the LB result is exact. However, for the ther-moelectric effects discussed in Sec. V which require thepresence of strong Coulomb interaction, an alternativemethod such as the ME approach is needed.We consider a single-level QD coupled to two leads (cid:96) ∈ {
A, B } (such as System 1 in Fig. 1 when tunnel- andCoulomb-decoupled from System 2). For simplicity, weassume energy-independent lead couplings γ (cid:96) ( (cid:15) ) = γ (cid:96) inthis case. The Hamiltonian of the QD reduces toˆ H dots = (cid:15) ˆ c † ˆ c , (13)with states labeled by the occupancy, | n (cid:105) ∈ {| (cid:105) , | (cid:105)} .In the LB formalism, the electric current and heat cur-rent going into lead A are given by, I LB A = − eh (cid:90) d(cid:15) T ( (cid:15) )[ f B ( (cid:15) ) − f A ( (cid:15) )] , (14)and J LB A = 1 h (cid:90) d(cid:15) ( (cid:15) − µ A ) T ( (cid:15) )[ f B ( (cid:15) ) − f A ( (cid:15) )] , (15)respectively. For a non-interacting single-level QD thetransmission function T ( (cid:15) ) is T ( (cid:15) ) = γ A γ B ( (cid:15) − (cid:15) ) + ( γ/ , (16)where γ = γ A + γ B and we have omitted the tunneling-induced energy shift which is not captured by the T -matrix approach.The transport currents calculated with the two ap-proaches with a finite bias and temperature difference( T B = 2 T A ≡ T ) between the leads are plotted inFigs. 2(a) and 2(b) as a function of the gate voltage fortwo different lead coupling strengths. To demonstrate theimportance of cotunneling processes, we have includedME results based on sequential tunneling only (black dot-ted curves) which do not depend on γ (cid:96) in the units shown,as well as sequential plus cotunneling (dashed curves). The results based purely on sequential tunneling differsignificantly from the LB results unless γ (cid:96) (cid:28) k B T . How-ever, for γ (cid:96) < k B T , the ME results with cotunnelingare in excellent agreement with the LB formalism. For γ (cid:96) > k B T which is outside the regime of validity of theME approach, the two approaches deviate, as expected.In the following discussion of thermoelectric effects, theheat current is of particular interest. As seen in Fig.2(b), when the dot level is above the electrochemical po-tential in lead A , the heat current becomes negative (forsufficiently small lead coupling strength). In this case,electrons above the electrochemical potential tunnel outof the lead and thereby cool the lead [cf. Eq. (10)].Such cooling mechanisms due to energy-selective tun-neling have been confirmed experimentally in metallicQD systems. The energy-selective tunneling gives riseto an asymmetry in the energy dissipation between thesource and drain leads which was recently observed inmolecular junctions. V. THERMOELECTRIC EFFECTS INCOULOMB-COUPLED QDS
In the remaining part of the paper, we study the ther-moelectric properties of the system illustrated in Fig. 1,i.e. two single-level QDs with QD1 tunnel-coupled toleads A and B and QD2 tunnel-coupled to lead C . TheCCQD system is described by the Hamiltonianˆ H dots = (cid:15) ˆ c † ˆ c + (cid:15) ˆ c † ˆ c + U ˆ n ˆ n , (17)where we have used the simplified notation U ≡ U , and the occupation states are | m (cid:105) = | n n (cid:105) ∈{| (cid:105) , | (cid:105) , | (cid:105) , | (cid:105)} . We consider situations where a J A [ γ ‘ ( ~ β ) − ] I A [ − e γ ‘ ~ − ] V [( eβ ) − ] (b)-100 5 10-5 V [( eβ ) − ]-10 MELBMELB (0 . . . . γ ‘ β )MethodME, Seq. FIG. 2. Comparison of the electric current (a) and heat cur-rent (b) calculated with the ME and LB approaches. Cur-rents are plotted as function of gate voltage V for two differ-ent lead coupling strengths γ A = γ B = γ (cid:96) (energy indepen-dent). The ME result including only sequential tunneling isshown for reference (black dotted), and the vertical dashedlines mark the alignment of the dot level with the electro-chemical potentials of lead A (left) and B (right). Parame-ters: T B = 2 T A ≡ / ( k B β ), µ A = 3 β − , µ B = − β − , and η = 10 − β − . source-drain bias V is applied to System 1, µ A = µ + eV / µ B = µ − eV / µ = 0 as reference).As pointed out above, we here allow for energy-dependent lead couplings. For bias voltages and tem-perature differences small compared to the energy scaleat which the lead couplings vary, it suffices to considerthe expansion of the lead couplings around their value at µ , γ (cid:96) ( (cid:15) ) = γ (cid:96) + ( (cid:15) − µ ) ∂γ (cid:96) , (18)where γ (cid:96) = γ (cid:96) ( µ ), ∂γ (cid:96) ≡ ∂γ (cid:96) ( (cid:15) ) ∂(cid:15) | (cid:15) = µ . A. Current and energy exchange
In Fig. 3(a) we show the electric current through QD1, I ≡ I A = − I B , at low temperature k B T (cid:96) = 10 − U (forillustrative convenience) and bias eV = 0 . U as a func-tion of gate detuning V − V and total gating V + V in the vicinity of the honeycomb vertex of the stabilitydiagram. Here, we initially assume energy-independentlead couplings which is sufficient to get an overall un-derstanding of the behavior of the system. The largecurrent near the degeneracy lines defined by ∆ , = 0and ∆ , = 0 is due to sequential tunneling processes.Away from these degeneracy lines where sequential tun-neling is exponentially suppressed, cotunneling processesgive rise to a weak background current. At the degen-eracy line ∆ , = 0 connecting the two triple pointsat V = V = 0 , U , respectively, nonlocal cotunnelingprocesses are responsible for the enhanced cotunnelingcurrent.The heat currents which accompany the electric cur-rent are shown in Figs. 3(b)–3(d) for different tempera-tures in the leads. Figure 3(b) shows the heat currentin lead A for k B T (cid:96) = 0 . U . Along the degeneracy lineswhere ∆ , = 0 and ∆ , = 0 and only the occu-pation of QD1 fluctuates, the heat current shows a be-havior similar to the one in Fig. 2(b) for a single-levelQD. However, at the center of the stability diagram,Coulomb-mediated energy exchange due to the strongCoulomb interaction between the QDs becomes signifi-cant. This manifests itself in a cooling of System 1 insidethe region bounded by the solid lines at the center ofFig. 3(b) (notice that the color scale is dominated bythe heat current with larger magnitude outside this re-gion). From the heat current in lead C shown in Fig. 3(c),the cooling of System 1 is seen to be at the cost ofheating System 2. This Coulomb-mediated energy ex-change between the two QD systems occurs in spite ofthe fact that no electrons are exchanged, and is the driv-ing force behind demon-induced cooling, energy har-vesting, and Coulomb drag. A simple analytical result for the energy exchange canbe found when considering sequential tunneling processesonly (indicated by the superscript s ). In this case, the to-tal heat currents in System 1, J s ≡ J sA + J sB , and System (c) J C [10 − U / ~ ] − V − V [ U/e ] V + V [ U / e ] −
11 6420 (b) J A [10 − U / ~ ] − V − V [ U/e ] V + V [ U / e ] −
11 0 JA =0 JB =0 (a) log[ I/ ( eγ A/B / (2 ~ ))] − V − V [ U/e ] V + V [ U / e ] − − ∆ , = ∆ , = ∆ , = − − (d) − − T ‘ [ U/k B ] J C [ U / ~ ] − −
40 10 − − − − − FIG. 3. Electric current and heat currents. (a) Electric cur-rent in System 1 as function of gate detuning V − V and totalgating V + V at low temperature, k B T (cid:96) = 10 − U . (b) Heatcurrent in lead A , J A , at high temperature, k B T (cid:96) = 10 − U (contours indicate where J A and J B are zero). (c) Heat cur-rent in lead C , J C , for k B T (cid:96) = 10 − U . (d) J C as functionof temperature with (solid) and without (dashed) cotunnelingfor the gate configurations marked in (c): eV , = 0 . U (blackcircle) and eV = 0 . U , eV = 0 . U (blue triangle). In plots(a)–(c), the degeneracy lines of the honeycomb vertex are in-dicated with dotted lines. Parameters: γ A/B ( (cid:15) ) = 10 − U , γ C ( (cid:15) ) = 10 − U , and eV = 0 . U . J s ≡ J sC , become J s = Uτ s (Γ + − Γ − ) + µ A − µ B e I s , (19a) J s = Uτ s (Γ − − Γ + ) , (19b)where Γ − ≡ Γ , Γ , Γ , Γ , , Γ + ≡ Γ , Γ , Γ , Γ , . The factor τ s depends onthe various sequential tunneling rates, however, ismerely a normalization factor and is not reproducedhere. The first two terms proportional to U in Eq. (19)describe the energy exchange, whereas the last term inEq. (19a) describes the contribution from Joule heatingin System 1. The direction of the energy transfer isdetermined by the sign of Γ − − Γ + . It is thereforeconvenient to consider the ratioΓ − Γ + = Ω e U ( β − β ) , (20)which describes whether energy is transferred from Sys-tem 1 to 2 (Γ − / Γ + >
1) or vice versa (Γ − / Γ + < On the right-hand side of (20), we have taken β A/B = β and β C = β , and expressed the ratio in terms of an ex-ponential factor, which depends on the temperature inSystem 1 and System 2, andΩ ≡ ( γ A f A + γ B f B )( γ A f A e − β µ A + γ B f B e − β µ B )( γ A f A + γ B f B )( γ A f A e − β µ A + γ B f B e − β µ B ) , (21)which depends on the temperature and bias in System1 only. The subscript 0 (1) in Eq. (21) indicates thatthe corresponding function is evaluated at ∆ , (∆ , )[see Eqs. (6)–(7)].The exponential factor in (20) shows that a tempera-ture gradient between the two QD systems can generate anet heat flow from the hot to the cold system. This is themechanism behind the heat engine studied in Ref. 11. Onthe other hand, a closer inspection of the Ω factor revealsthat it is, in fact, possible to generate a net heat flow inthe opposite direction, i.e. from the cold to the hot sys-tem, and this is the cause of the negative heat currentat the center of Fig. 3(b). This so-called demon-inducedcooling effect will be discussed further in Sec. V B below.When the applied bias and temperature are small com-pared to the inter-dot Coulomb interaction, eV, k B T (cid:28) U , cotunneling processes start to dominate the heat cur-rents. This is demonstrated in Fig. 3(d) which showsthe heat current J C as a function of temperature forthe two different gate tunings marked with symbols inFig. 3(c). Considering sequential tunneling only (dashedcurves), the heat current is quenched at k B T (cid:28) U as Γ , and Γ , in Γ − become exponentially sup-pressed. This can also be understood from the illustra-tion in Fig. 4(a) which shows the sequence of sequen-tial tunneling processes corresponding to Γ − . However,nonlocal cotunneling processes allow the system to fluc-tuate between the two states 10 ↔
01, as illustratedin Fig. 4(b), and thereby transfer heat between thesystems. The nonlocal cotunneling channel is open for | ∆ , | (cid:46) max {| eV / | , k B T } , and the associated heatcurrent is thus also suppressed at low temperature when∆ , (cid:54) = 0 as illustrated by the blue curve (triangle) inFig. 3(d). For zero detuning ∆ , = 0 (circle), the non-local cotunneling rates, and hence also the heat current,saturate at k B T (cid:28) eV . In Sec. V B, we demonstratethat nonlocal cotunneling processes have a significant ef-fect on the demon-induced cooling mechanism. B. Demon-induced cooling
The effect of cooling System 1 at the cost of heat-ing System 2 has recently been discussed in context ofa Maxwell’s demon where System 2 plays the role ofthe demon which performs the necessary feedback tocool System 1.
To maximize the achievable coolingpower for refrigeration purposes, large tunneling rates, γ (cid:96) ( (cid:15) ) ∼ k B T, U , are essential [cf. Eq. (19)]. How-ever, large tunneling rates increase the contribution fromhigher-order tunneling processes, thus emphasizing theimportance of including cotunneling processes in quanti-tative analyses of the cooling power. U Γ −→ A ←− C ,
10 Γ −→ C ←− B , U Γ ←− C , ←− B , −→ C , − J [ − U / ~ ] J s J c J ( a ) ( b ) . . V [ U/e ] ( c ) Γ −→ A , . . . . FIG. 4. Cooling cycle and effect of cotunneling. (a) Sequenceof sequential tunneling processes which cools System 1. Thepositions of the dot levels when the other dot is empty (occu-pied) is illustrated with solid (dotted) lines. (b) Sequenceof nonlocal cotunneling processes. (c) Heat current J asfunction of bias voltage. The individual contributions fromsequential ( J s ) and cotunneling ( J c ) are also shown. Param-eters: eV = eV = U/ γ A/B ( (cid:15) ) = 10 − U , γ C ( (cid:15) ) = 10 − U ,and k B T = 0 . U . In the following, we consider the case of uniform tem-perature T (cid:96) ≡ T whereby the exponential factor in (20)becomes unity. This allows us to focus on the Ω factorin the optimization of the performance. Equation (19)shows that the cooling mechanism is governed by Γ − since, as illustrated in Fig. 4(a), in a full sequential cyclean amount of energy U is transferred from System 1 toSystem 2 thereby cooling System 1. In the following, wediscuss how to increase the cooling power by maximiz-ing the success rate for completing the cooling cycle inFig. 4(a).
1. Cotunneling limitations
Although the cycle of nonlocal cotunneling processesillustrated in Fig. 4(b) gives the same net transfer ofelectrons as the sequential tunneling cycle in Fig. 4(a),the net energy transfer is different for the two cases. Asillustrated, in a cotunneling process also electrons below(above) the electrochemical potential can tunnel out oflead A (into lead B ), and thus reduce the demon-inducedcooling effect.In Fig. 4(c), we show the heat current J together withits individual contributions from sequential ( J s ) and co-tunneling ( J c ) processes. Overall, System 1 cools atlow bias, while at higher bias, Joule heating becomesdominant. The minimum in J as a function of biasvoltage is referred to as the maximum cooling power, J , max ≡ min J ( V ). As the figure shows, cotunnelingreduces the maximum cooling power.Figure 5 shows how the maximum cooling power J , max scales with the lead coupling strengths. As the figuredemonstrates, the rates must satisfy γ C > γ A/B to en- -4-10-6-8 log( γ A/B /U ) − − − − − − l og ( γ C / U ) (b) log[ − J , max / ( U / ~ )] − log( γ A/B /U ) − − − − − − − l og ( γ C / U ) (a) log[ − J s , max / ( U / ~ )] -4-10 − − − − − − − − − -6-8 FIG. 5. Maximum cooling power, J , max , as function of thelead coupling strengths for energy-independent couplings. (a)Sequential tunneling, and (b) sequential plus cotunneling. Pa-rameters: eV = eV = U/ k B T = 0 . U . sure that System 2 acts sufficiently fast to perform thedesired feedback such that the cooling cycle in Fig. 4(a) iscompleted when an electron tunnels between lead A and B . In the region of large cooling power, cotunnelingprocesses start to become important, and hence there isa trade-off between sequential tunneling which improvesthe cooling effect, and nonlocal cotunneling which lim-its the effect. In addition, the area in the lead couplingparameter space where refrigeration is possible is also re-duced when cotunneling is included.
2. Performance boosting
Here we demonstrate that energy-dependent lead cou-plings can enhance the demon-induced cooling power sig-nificantly. We restrict the discussion to lead couplingswith a linear energy dependence [cf. Eq. (18)].By inspecting the Ω factor in Eq. (21), we find thatfor µ A > µ B , the configuration illustrated in the inset ofFig. 6 where γ A , γ B are reduced compared to γ A , γ B ,boosts the Ω factor (and thereby Γ − / Γ + ). This resultsin an enhancement of the cooling power by suppressingdirect tunneling between lead A and B via two sequentialtunneling processes which contributes to Joule heating,while at the same time promoting the processes of thecooling cycle in Fig. 4(a).In Fig. 6 we show the maximum cooling power as afunction of temperature for different situations for theenergy dependence of the lead couplings, from the top(black) curve showing the result for energy independent lead couplings, to increasing energy dependence, i.e. in-creasing | ∂γ A/B | , towards the bottom (light blue) curve.When tuning the energy dependence of the lead cou-plings, a significant enhancement of the cooling poweris achieved. Again, the effect of cotunneling processes isto reduce the attainable cooling power (solid lines) rela-tive to the cooling power obtained when only consideringsequential tunneling processes (dashed lines). VI. CONCLUSIONS
In summary, we have studied thermoelectric effects inCCQD systems with a T -matrix based master-equationapproach for the calculation of charge and heat currents.Importantly, our method (i) treats incoherent sequen-tial tunneling processes and coherent cotunneling pro-cesses on equal footing, and (ii) can account for energy-dependent tunnel couplings to the leads. Both are essen-tial for quantitative predictions and optimization of thethermoelectric properties of CCQDs.To benchmark the master-equation method, we consid-ered a non-interacting single-level QD coupled to sourceand drain leads for which the Landauer-B¨uttiker formal-ism is exact. In the regime of validity of our method, i.e.small tunnel couplings to the leads, γ < k B T , we demon-strated excellent agreement with the results from theLandauer-B¨uttiker method when cotunneling processesare included in the master equation.Furthermore, we studied the effect of cotunneling pro-cesses and energy-dependent lead couplings on the ther-moelectric properties of a CCQD system consisting of twoQDs exhibiting a Maxwell’s demon-like cooling mecha-nism. First of all, we showed that cotunneling pro-cesses reduce the cooling effect since cotunneling pro-cesses do not share the delicate energy selectivity in-herent to sequential tunneling processes. This results ina significant reduction of the achievable cooling powercompared to the sequential tunneling result when thelead couplings are increased to maximize the coolingpower from sequential tunneling processes. Secondly,we demonstrated that it is possible to boost the cool-
A Bk B T [ U ] J , m a x [ − U / ~ ] − − −
10 0 . . . J , max J s , max . . . FIG. 6. Performance boosting with energy-dependent leadcouplings. Maximum cooling power as function of tempera-ture for different lead coupling strengths: ∂γ A = − ∂γ B = xγ A/B /U (sketched in the inset), with x = 0 (black) to x = 1(light blue) in steps of 0 .
2. The full (dashed) lines showthe result obtained with (without) cotunneling. Parameters: γ C ( (cid:15) ) = 10 − U , γ A/B = 10 − U , eV = eV = U/
2, and η = 10 − U . ing power significantly via other means by introduc-ing energy-dependent lead couplings and properly tun-ing their energy dependence. In this case, we showedthat cotunneling still reduces the cooling power signifi-cantly, thus emphasizing the importance of cotunnelingprocesses in quantitative analyses.Applying the methodology to other mesoscopic sys-tems allows for testing of new thermoelectric device ideasbeyond sequential tunneling estimates, as well as for im-proved comparison with experiments. ACKNOWLEDGMENTS
We would like to thank J. P. Pekola, M. Leijnse,C. Timm, and N. M. Gergs for valuable discussions.K.K. acknowledges support from the European Union’sHorizon 2020 research and innovation programme underthe Marie Sklodowska-Curie grant agreement no. 713683(COFUNDfellowsDTU). The Center for NanostructuredGraphene (CNG) is sponsored by the Danish ResearchFoundation, Project DNRF103.
Appendix A: Cotunneling rates and regularizationprocedure
The rate for elastic cotunneling through a single-levelQD is given by˜Γ −→ (cid:96) ←− (cid:96) (cid:48) mm = (cid:90) d(cid:15) π (cid:126) γ (cid:96) ( (cid:15) ) γ (cid:96) (cid:48) ( (cid:15) ) f (cid:96) ( (cid:15) ) ¯ f (cid:96) (cid:48) ( (cid:15) ) (cid:12)(cid:12)(cid:12)(cid:12) vm ± (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (A1)where v refers to the virtually occupied intermediatestate created in the process where an initially empty levelis filled (+ (cid:15) ) or an initially filled level is emptied ( − (cid:15) ).In pair-cotunneling processes, two electrons tunnel si-multaneously out of (into) the QD system and into (outof) the leads (cid:96) and (cid:96) (cid:48) . The rate for such processes takesthe form˜Γ ←− (cid:96) ←− (cid:96) (cid:48) mn = (cid:90) d(cid:15) π (cid:126) γ (cid:96) ( (cid:15) ) γ (cid:96) (cid:48) (∆ nm − (cid:15) ) ¯ f (cid:96) ( (cid:15) ) ¯ f (cid:96) (cid:48) (∆ nm − (cid:15) ) × (cid:12)(cid:12)(cid:12)(cid:12) vm − (cid:15) + iη + 1∆ v (cid:48) n + (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (A2)where v ( v (cid:48) ) refers to the virtually occupied intermediatestate in a process where an electron initially tunnels fromthe QD system and into lead (cid:96) ( (cid:96) (cid:48) ). Similarly,˜Γ −→ (cid:96) −→ (cid:96) (cid:48) mn = (cid:90) d(cid:15) π (cid:126) γ (cid:96) ( (cid:15) ) γ (cid:96) (cid:48) (∆ mn − (cid:15) ) f (cid:96) ( (cid:15) ) f (cid:96) (cid:48) (∆ mn − (cid:15) ) × (cid:12)(cid:12)(cid:12)(cid:12) vn − (cid:15) + iη + 1∆ v (cid:48) m + (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (A3) where v ( v (cid:48) ) refer to the virtually occupied intermediatestate in a process where an electron initially tunnels fromlead (cid:96) (cid:48) ( (cid:96) ) and into the QD system.The bare cotunneling rates are formally divergent inthe limit η →
0. The divergence stems from factors in-volving | x + iη | − , x, η ∈ R . Using that (cid:12)(cid:12)(cid:12)(cid:12) x + iη (cid:12)(cid:12)(cid:12)(cid:12) → πη δ ( x ) + P x , η → + , (A4)where P denotes the principle value, we can identify thedivergent contributions, e.g. from Eq. (8)˜Γ −→ (cid:96) ←− (cid:96) (cid:48) mn → (cid:126) η (cid:16) Γ −→ (cid:96) mv Γ ←− (cid:96) (cid:48) vn + Γ ←− (cid:96) (cid:48) mv (cid:48) Γ −→ (cid:96) v (cid:48) n (cid:17) + Γ −→ (cid:96) ←− (cid:96) (cid:48) mn , (A5)where Γ −→ (cid:96) ←− (cid:96) (cid:48) mn denotes the regularized cotunneling rate, andwe have used that the cross-terms from the absolutesquared in Eq. (8) do not contribute to any divergences.The divergent contribution is proportional to products oftwo sequential tunneling rates. These correspond to twoenergy-conserving (sequential) transitions which can beidentified with the intermediate processes in the cotun-neling process. The sum is over the possible sequences ofintermediate transitions. Similarly, for the cotunnelingheat rates, e.g. Eq. (12)˜ W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96),mn → (cid:126) η (cid:104) W −→ (cid:96) (cid:96),mv Γ ←− (cid:96) (cid:48) vn + Γ ←− (cid:96) (cid:48) mv (cid:48) W −→ (cid:96) (cid:96),v (cid:48) n (cid:105) + W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96),mn , (A6)or the corresponding heat rate in lead (cid:96) (cid:48) ˜ W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96) (cid:48) ,mn → (cid:126) η (cid:104) Γ −→ (cid:96) mv W ←− (cid:96) (cid:48) (cid:96) (cid:48) ,vn + W ←− (cid:96) (cid:48) (cid:96) (cid:48) ,mv (cid:48) Γ −→ (cid:96) v (cid:48) n (cid:105) + W −→ (cid:96) ←− (cid:96) (cid:48) (cid:96) (cid:48) ,mn . (A7)We apply the regularization scheme in Ref. 29 and sub-tract these terms scaling as η − .In the case of identical temperatures in the leads, usingthe identity f ( (cid:15) )[1 − f ( (cid:15) )] = n ( (cid:15) − (cid:15) )[ f ( (cid:15) ) − f ( (cid:15) )],where f ( (cid:15) ) is the Fermi-Dirac distribution and n ( (cid:15) ) is theBose-Einstein distribution, the cotunneling rates can bewritten in the form I = (cid:90) ∞−∞ d(cid:15)P ( (cid:15) ) (cid:104) f (cid:96) (cid:48) ( (cid:15) ) − f (cid:96) ( (cid:15) + ∆ ) (cid:105) × (cid:12)(cid:12)(cid:12)(cid:12) k (cid:15) − ∆ + iη + k ∆ − (cid:15) + iη (cid:12)(cid:12)(cid:12)(cid:12) , (A8)where P ( (cid:15) ) is assumed to be a polynomial, P ( (cid:15) ) = (cid:80) ni =0 c i (cid:15) n , of maximum order n = 2 for k − k (cid:54) = 0and n = 4 for k − k = 0 to ensure that the result belowis well-defined. The derivation is in line with the one inRef. 27, and the integral becomes I = k P (cid:48) (∆ )Re (cid:2) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:3) + k β π P (∆ )Im (cid:104) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:105) + k P (cid:48) (∆ )Re (cid:2) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:3) + k β π P (∆ )Im (cid:104) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:105) − k k ∆ − ∆ (cid:0) P (∆ )Re (cid:2) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:3) − P (∆ )Re (cid:2) ψ − (cid:96) (cid:48) (∆ ) − ψ − (cid:96) (∆ + ∆ ) (cid:3)(cid:1) + R + O ( η − ) + O ( η ) , (A9) where ψ (1) ± (cid:96) ( (cid:15) ) ≡ ψ (1) (cid:18) ± i β π ( (cid:15) − µ (cid:96) ) (cid:19) , (A10)with ψ ( ψ ) being the digamma (trigamma) function, and R = (cid:26) c ( µ (cid:96) (cid:48) − µ (cid:96) + ∆ )( k − k ) , k − k (cid:54) = 0 ,c ( µ (cid:96) (cid:48) − µ (cid:96) + ∆ ) k (∆ − ∆ ) , k − k = 0 . (A11) The term O ( η − ) is omitted by regularization before let-ting η →
0. For k B T < γ (outside the regime of valid-ity), the failure of the approach is seen as a logarithmic divergence of the digamma functions near the degeneracypoints.In studies of thermoelectric effects where different leadtemperatures as well as more general energy-dependenceof the lead couplings become relevant, one must turn toa numerical procedure. In this case, we evaluate the co-tunneling integrals numerically with a small but finite η ,and subsequently subtract contributions of order η − asshown in e.g. Eqs. (A5)–(A7). In particular, we haveapplied the numerical procedure in Fig. 2 and Fig. 6,and stated the values of η in the figure caption. ∗ [email protected] † [email protected] J. P. Pekola, O.-P. Saira, V. F. Maisi, A. Kemppinen,M. M¨ott¨onen, Y. A. Pashkin, and D. V. Averin, “Single-electron current sources: Toward a refined definition of theampere,” Rev. Mod. Phys. , 1421–1472 (2013). F. Giazotto, T. T. Heikkil¨a, A. Luukanen, A. M. Savin,and J. P. Pekola, “Opportunities for mesoscopics in ther-mometry and refrigeration: Physics and applications,”Rev. Mod. Phys. (2006). M. S. Dresselhaus, G. Chen, M. Y. Tang, R. G. Yang,H. Lee, D. Z. Wang, Z. F. Ren, J.-P. Fleurial, andP. Gogna, “New directions for low-dimensional thermoelec-tric materials,” Adv. Mater. , 1043–1053 (2007). H. Thierschmann, R. S´anchez, B. Sothmann, H. Buhmann,and L. W. Molenkamp, “Thermoelectrics with Coulomb-coupled quantum dots,” C. R. Physique , 1109 – 1122(2016). D. Bischoff, M. Eich, O. Zilberberg, C. R¨ossler, T. Ihn,and K. Ensslin, “Measurement back-action in stackedgraphene quantum dots,” Nano. Lett. , 6003 (2015). A. J. Keller, J. S Lim, D. S´anchez, R. L´opez, S. Amasha,J. A. Katine, H. Shtrikman, and D. Goldhaber-Gordon,“Cotunneling drag effect in Coulomb-coupled quantumdots,” Phys. Rev. Lett. , 066602 (2016). A. Hamo, A. Benyamini, I. Shapir, I. Khivrich, J. Waiss-man, K. Kaasbjerg, Y. Oreg, F. von Oppen, and S. Ilani,“Electron attraction mediated by Coulomb repulsion,” Na-ture , 395 (2016). H. Thierschmann, R. S´anchez, B. Sothmann, F. Arnold,C. Heyn, W. Hansen, H. Buhmann, and L. W.Molenkamp, “Three-terminal energy harvester with cou-pled quantum dots,” Nature Nanotech. , 854 (2015). J. V. Koski, A. Kutvonen, I. M. Khaymovich, T. Ala-Nissila, and J. P. Pekola, “On-chip Maxwell’s demon as aninformation-powered refrigerator,” Phys. Rev. Lett. ,260602 (2015). P. Strasberg, G. Schaller, T. Brandes, and M. Espos- ito, “Thermodynamics of a physical model implementinga Maxwell demon,” Phys. Rev. Lett. , 040601 (2013). R. S´anchez and M. B¨uttiker, “Optimal energy quanta tocurrent conversion,” Phys. Rev. B , 085428 (2011). B. Sothmann, R. S´anchez, and A. N. Jordan, “Thermo-electric energy harvesting with quantum dots,” Nanotech-nology , 032001 (2014). R. S´anchez, H. Thierschmann, and L. W. Molenkamp,“All-thermal transistor based on stochastic switching,”(2017), arXiv:1701.00382. J. V. Koski and J. P. Pekola, “Maxwell’s demons realized inelectronic circuits,” C. R. Physique , 1130 – 1138 (2016). G. Benenti, G. Casati, K. Saito, and R. S. Whitney, “Fun-damental aspects of steady-state conversion of heat to workat the nanoscale,” (2016), arXiv:1608.05595. F. Haupt, M. Leijnse, H. L. Calvo, L. Classen,J. Splettstoesser, and M. R. Wegewijs, “Heat, molecularvibrations, and adiabatic driving in non-equilibrium trans-port through interacting quantum dots,” Phys. Status So-lidi B , 2315–2329 (2013). M. Leijnse, M. R. Wegewijs, and K. Flensberg, “Nonlin-ear thermoelectric properties of molecular junctions withvibrational coupling,” Phys. Rev. B , 045412 (2010). J. Arg¨uello-Luengo, D. S´anchez, and R. L´opez, “Heatasymmetries in nanoscale conductors: The role of deco-herence and inelasticity,” Phys. Rev. B , 165431 (2015). D. S´anchez and R. L´opez, “Nonlinear phenomena in quan-tum thermoelectrics and heat,” C. R. Physique , 1060 –1071 (2016). N. M. Gergs, C. B. M. H¨orig, M. R. Wegewijs, andD. Schuricht, “Charge fluctuations in nonlinear heat trans-port,” Phys. Rev. B , 201107(R) (2015). Kevin Marc Seja, Gediminas Kirˇsanskas, Carsten Timm,and Andreas Wacker, “Violation of Onsager’s theorem inapproximate master equation approaches,” Phys. Rev. B , 165435 (2016). A.-M. Dar´e and P. Lombardo, “Powerful Coulomb dragthermoelectric engine,” (2017), arXiv:1704.04064. H. Bruus and K. Flensberg,
Many-body Quantum Theoryin Condensed Matter Physics (Oxford University Press,2004). J. Waissman, M. Honig, S. Pecker, A. Benyamini,A. Hamo, and S. Ilani, “Realization of pristine and lo-cally tunable one-dimensional electron systems in carbonnanotubes,” Nature Nanotech. , 569 (2013). Y. Zhang, G. Lin, and J. Chen, “Three-terminal quantum-dot refrigerators,” Phys. Rev. E , 052118 (2015). R. S´anchez, R. L´opez, D. S´anchez, and M. B¨uttiker,“Mesoscopic Coulomb drag, broken detailed balance,and fluctuation relations,” Phys. Rev. Lett. , 076801(2010). K. Kaasbjerg and A.-P. Jauho, “Correlated Coulomb dragin capacitively coupled quantum-dot structures,” Phys.Rev. Lett. , 196801 (2016). J. S. Lim, R. L´opez, and D. S´anchez, “Engineering dragcurrents in Coulomb coupled quantum dots,” (2016),arXiv:1612.06627. M. Turek and K. A. Matveev, “Cotunneling thermopowerof single electron transistors,” Phys. Rev. B , 115332(2002). J. Koch, F. von Oppen, Y. Oreg, and E. Sela, “Ther-mopower of single-molecule devices,” Phys. Rev. B ,195107 (2004). C. Timm, “Tunneling through molecules and quantumdots: Master-equation approaches,” Phys. Rev. B ,195416 (2008). J. K¨onig, J. Schmid, H. Schoeller, and G. Sch¨on, “Reso-nant tunneling through ultrasmall quantum dots: Zero-bias anomalies, magnetic-field dependence, and boson-assisted transport,” Phys. Rev. B , 16820–16837 (1996). A. Thielmann, M. H. Hettler, J. K¨onig, and G. Sch¨on,“Cotunneling current and shot noise in quantum dots,”Phys. Rev. Lett. , 146806 (2005). J. N. Pedersen and A. Wacker, “Tunneling throughnanosystems: Combining broadening with many-particlestates,” Phys. Rev. B , 195330 (2005). C. Timm, “Time-convolutionless master equation for quan-tum dots: Perturbative expansion to arbitrary order,”Phys. Rev. B , 115416 (2011). S. Amasha, A. J. Keller, I. G. Rau, A. Carmi, J. A. Ka-tine, H. Shtrikman, Y. Oreg, and D. Goldhaber-Gordon,“Pseudospin-resolved transport spectroscopy of the Kondoeffect in a double quantum dot,” Phys. Rev. Lett. ,046604 (2013). J. Koch, M. E. Raikh, and F. von Oppen, “Pair tunnel-ing through single molecules,” Phys. Rev. Lett. , 056803(2006). M. Leijnse, M. R. Wegewijs, and M. H. Hettler, “Pair tun-neling resonance in the single-electron transport regime,”Phys. Rev. Lett. , 156803 (2009). S. Koller, M. Grifoni, M. Leijnse, and M. R. Wegewijs,“Density-operator approaches to transport through inter-acting quantum dots: Simplifications in fourth order per-turbation theory,” Phys. Rev. B , 235307 (2010). Notice that for the particular pair-cotunneling processeswith (cid:96) = (cid:96) (cid:48) one should include a factor of two in the cur-rent, however, for single-level QDs discussed here such pro-cesses do not contribute. M. F. Ludovico, M. Moskalets, D. S´anchez, and L. Ar-rachea, “Dynamics of energy transport and entropy pro-duction in ac-driven quantum electron systems,” Phys.Rev. B , 035436 (2016). J. Koch, F. von Oppen, and A. V. Andreev, “Theory of thefranck-condon blockade regime,” Phys. Rev. B , 205438(2006). H. Haug and A.-P. Jauho,
Quantum Kinetics in Transportand Optics of Semiconductors (Springer Berlin Heidelberg,2008). A. V. Feshchenko, J. V. Koski, and J. P. Pekola, “Exper-imental realization of a Coulomb blockade refrigerator,”Phys. Rev. B , 201407 (2014). W. Lee, K. Kim, W. Jeong, L. A. Zotti, F. Pauly, J. C.Cuevas, and P. Reddy, “Heat dissipation in atomic-scalejunctions,” Nature , 209–U103 (2013). To ensure that the lead coupling strengths are positive, alinear expansion is only appropriate when the bias windowor the thermal window exponentially suppress the contri-bution to the cotunneling integrals at energies where thelead coupling becomes negative. In the numerical calcula-tion we take the absolute value of the lead couplings. W. G. van der Wiel, S. De Franceschi, J. M. Elzerman,T. Fujisawa, S. Tarucha, and L. P. Kouwenhoven, “Elec-tron transport through double quantum dots,” Rev. Mod.Phys. , 1–22 (2002). R. S´anchez and M. B¨uttiker, “Detection of single-electronheat transfer statistics,” EPL , 47008 (2012). A. Kutvonen, J. Koski, and T. Ala-Nissila, “Thermody-namics and efficiency of an autonomous on-chip Maxwell’sdemon,” Scientific Reports6