Thermoelectrics in Coulomb-coupled quantum dots: Cotunneling and energy-dependent lead couplings

TL;DR

This paper develops a master-equation approach to study thermoelectric effects in Coulomb-coupled quantum dots, highlighting the significance of cotunneling and energy-dependent lead couplings in optimizing heat and cooling power.

Contribution

It introduces a perturbative T-matrix master-equation method that accurately includes cotunneling and energy-dependent lead couplings, advancing the modeling of thermoelectric transport in quantum-dot systems.

Findings

Cotunneling dominates heat currents at low temperature and bias.

Energy-dependent lead couplings can significantly enhance cooling power.

Higher-order processes reduce the effectiveness of Maxwell's demon cooling.

Abstract

We study thermoelectric effects in Coulomb-coupled quantum-dot (CCQD) systems beyond lowest-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 the charge and heat tunneling rates and transport currents. Applying the method to transport through a non-interacting single-level QD, we demonstrate excellent agreement with the Landauer-B{\"u}ttiker theory when higher-order (cotunneling) processes are included in the ME. Next, we study the effect of cotunneling and energy-dependent lead couplings on the heat currents in a system of two Coulomb-coupled QDs. Overall, we find that cotunneling processes (i) dominate the heat currents at low temperature and bias, and (ii) give rise to a pronounced reduction of the cooling power achievable with the recently demonstrated Maxwell's demon cooling mechanism. Furthermore, we demonstrate that the cooling power can be boosted significantly by carefully engineering the energy dependence of the lead couplings to filter out undesired transport processes. Our findings emphasize the importance of considering higher-order cotunneling processes as well as the advantage of engineered energy-dependent lead couplings in the optimization of the thermoelectric performance of Coulomb-coupled QD systems.

Equations

$$\hat{H}={\hat{H}}_{\text{dots}}+{\hat{H}}_{\text{leads}}+{\hat{H}}_T,$$

$${\hat{H}}_{\text{dots}}=\sum _{\delta }{ϵ}_{\delta }^{}{\hat{c}}_{\delta }^{†}{\hat{c}}_{\delta }^{}+\sum _{⟨\delta ,{\delta }^{\prime }⟩}{U}_{\delta {\delta }^{\prime }}^{}{\hat{n}}_{\delta }^{}{\hat{n}}_{{\delta }^{\prime }}^{},$$

$${\dot{p}}_m=\sum _{n\ne m}\left({\Gamma }_{nm}{p}_n-{\Gamma }_{mn}{p}_m\right),\sum _m{p}_m=1,$$

$${\tilde{\Gamma }}_{mn}=\frac{2\pi }{\hslash }\sum _{ij}|⟨j|⟨n|T|m⟩|i⟩{|}^2{\rho }_i\delta ({\Delta }_{mn}+E_{\text{leads},j}-E_{\text{leads},i}),$$

$$\hat{T}={\hat{H}}_T+{\hat{H}}_T\frac{1}{E_{\text{initial}}-{\hat{H}}_0+i\eta }\hat{T},$$

$${\Gamma }_{mn}^{\overrightarrow{\ell }}$$

$${\Gamma }_{mn}^{\overleftarrow{\ell }}$$

$${\tilde{\Gamma }}_{mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}$$

$$\times {\left|\frac{1}{{\Delta }_{vm}+ϵ+i\eta }+\frac{1}{{\Delta }_{v^{\prime }n}-ϵ+i\eta }\right|}^2,$$

$$I_{\ell }\equiv -e⟨\sum _k\frac{d{\hat{n}}_{\ell k}}{dt}⟩=-e\sum _{mn}{p}_m^{}\left({\Gamma }_{mn}^{\ell \leftarrow }-{\Gamma }_{mn}^{\ell \to }\right),$$

$$J_{\ell }\equiv ⟨\sum _k({ϵ}_{\ell k}-{\mu }_{\ell })\frac{d{\hat{n}}_{\ell k}}{dt}⟩=\sum _{mn}{p}_m^{}\left(W_{mn}^{\ell \leftarrow }-W_{mn}^{\ell \to }\right),$$

$$\begin{array}{rl}{W}_{\ell ,mn}^{\overrightarrow{\ell }}& =({\Delta }_{mn}-{\mu }_{\ell }){\Gamma }_{mn}^{\overrightarrow{\ell }},\\ W_{\ell ,mn}^{\overleftarrow{\ell }}& =({\Delta }_{nm}-{\mu }_{\ell }){\Gamma }_{mn}^{\overleftarrow{\ell }},\end{array}$$

$${\tilde{W}}_{\ell ,mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}$$

$$\times (ϵ-{\mu }_{\ell }){\left|\frac{1}{{\Delta }_{vm}+ϵ+i\eta }+\frac{1}{{\Delta }_{v^{\prime }n}-ϵ+i\eta }\right|}^2,$$

$${\hat{H}}_{\text{dots}}={ϵ}_1^{}{\hat{c}}_1^{†}{\hat{c}}_1^{},$$

$$I_A^{\text{LB}}=\frac{-e}{h}\int dϵT(ϵ)[f^B(ϵ)-f^A(ϵ)],$$

$$J_A^{\text{LB}}=\frac{1}{h}\int dϵ(ϵ-{\mu }_A)T(ϵ)[f^B(ϵ)-f^A(ϵ)],$$

$$T(ϵ)=\frac{{\gamma }^A{\gamma }^B}{(ϵ-{ϵ}_1{)}^2+(\gamma /2{)}^2},$$

$${\hat{H}}_{\text{dots}}={ϵ}_1^{}{\hat{c}}_1^{†}{\hat{c}}_1^{}+{ϵ}_2^{}{\hat{c}}_2^{†}{\hat{c}}_2^{}+U{\hat{n}}_1^{}{\hat{n}}_2^{},$$

$${\gamma }^{\ell }(ϵ)={\gamma }_0^{\ell }+(ϵ-{\mu }_0)\partial {\gamma }^{\ell },$$

$$J_1^s$$

$$J_2^s$$

$$\frac{{\Gamma }^{-}}{{\Gamma }^{+}}=\Omega e^{U({\beta }_2-{\beta }_1)},$$

$$\Omega \equiv \frac{({\gamma }_1^A{f}_1^A+{\gamma }_1^B{f}_1^B)({\gamma }_0^A{f}_0^A{e}^{-{\beta }_1{\mu }_A}+{\gamma }_0^B{f}_0^B{e}^{-{\beta }_1{\mu }_B})}{({\gamma }_0^A{f}_0^A+{\gamma }_0^B{f}_0^B)({\gamma }_1^A{f}_1^A{e}^{-{\beta }_1{\mu }_A}+{\gamma }_1^B{f}_1^B{e}^{-{\beta }_1{\mu }_B})},$$

$${\tilde{\Gamma }}_{mm}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}$$

$${\tilde{\Gamma }}_{mn}^{\overleftarrow{\ell }\overleftarrow{{\ell }^{\prime }}}$$

$$\times {\left|\frac{1}{{\Delta }_{vm}-ϵ+i\eta }+\frac{1}{{\Delta }_{v^{\prime }n}+ϵ+i\eta }\right|}^2,$$

$${\tilde{\Gamma }}_{mn}^{\overrightarrow{\ell }\overrightarrow{{\ell }^{\prime }}}$$

$$\times {\left|\frac{1}{{\Delta }_{vn}-ϵ+i\eta }+\frac{1}{{\Delta }_{v^{\prime }m}+ϵ+i\eta }\right|}^2,$$

$${\left|\frac{1}{x+i\eta }\right|}^2\to \frac{\pi }{\eta }\delta (x)+\mathcal{P}\frac{1}{x^2},\eta \to 0^{+},$$

$${\tilde{\Gamma }}_{mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}\to \frac{\hslash }{2\eta }\left({\Gamma }_{mv}^{\overrightarrow{\ell }}{\Gamma }_{vn}^{\overleftarrow{{\ell }^{\prime }}}+{\Gamma }_{m{v}^{\prime }}^{\overleftarrow{{\ell }^{\prime }}}{\Gamma }_{v^{\prime }n}^{\overrightarrow{\ell }}\right)+{\Gamma }_{mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}},$$

$${\tilde{W}}_{\ell ,mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}\to \frac{\hslash }{2\eta }\left[W_{\ell ,mv}^{\overrightarrow{\ell }}{\Gamma }_{vn}^{\overleftarrow{{\ell }^{\prime }}}+{\Gamma }_{m{v}^{\prime }}^{\overleftarrow{{\ell }^{\prime }}}{W}_{\ell ,v^{\prime }n}^{\overrightarrow{\ell }}\right]+W_{\ell ,mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}},$$

$${\tilde{W}}_{{\ell }^{\prime },mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}\to \frac{\hslash }{2\eta }\left[{\Gamma }_{mv}^{\overrightarrow{\ell }}{W}_{{\ell }^{\prime },vn}^{\overleftarrow{{\ell }^{\prime }}}+W_{{\ell }^{\prime },m{v}^{\prime }}^{\overleftarrow{{\ell }^{\prime }}}{\Gamma }_{v^{\prime }n}^{\overrightarrow{\ell }}\right]+W_{{\ell }^{\prime },mn}^{\overrightarrow{\ell }\overleftarrow{{\ell }^{\prime }}}.$$

Full text

I Introduction

The experimental progress in control of single-electron transport Pekola et al. (2013) has spurred interest in nanosystems which utilize the associated heat currents for thermoelectric applications Giazotto et al. (2006); Dresselhaus et al. (2007); Thierschmann et al. (2016). In particular, experiments with Coulomb-coupled quantum-dot (CCQD) systems have demonstrated a plethora of novel phenomena ranging from Coulomb drag Bischoff et al. (2015); Keller et al. (2016) and electron pairing Hamo et al. (2016) to extraordinary thermoelectric effects Thierschmann et al. (2015); Koski et al. (2015). This includes the realization of an energy harvester which converts a thermal gradient in a CCQD system into an electric current Thierschmann et al. (2015), as well as an autonomous Maxwell’s demon capable of cooling a current-carrying QD system at the cost of heating a “demon” QD system Koski et al. (2015).

In addition to the above, theoretical studies have predicted a wide range of novel thermoelectric effects in CCQD systems Strasberg et al. (2013); Sánchez and Büttiker (2011); Sothmann et al. (2014); Sánchez et al. (2017). The mechanisms behind these effects rely on the presence of a strong Coulomb interaction 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 thermoelectric properties of CCQD systems Thierschmann et al. (2016); Koski and Pekola (2016) and opens the opportunity to test fundamental thermodynamic aspects of heat transport in interacting nanoscale systems driven out of equilibrium Benenti et al. (2016).

While the operation principles of the above-mentioned effects are governed by incoherent electron tunneling (sequential tunneling) processes between the leads and the QDs Strasberg et al. (2013); Sánchez and Büttiker (2011); Sothmann et al. (2014); Thierschmann et al. (2015); Koski et al. (2015); Sánchez et al. (2017), the importance of coherent higher-order tunneling (cotunneling) processes for the nonlinear heat transport remains largely unexplored Koski et al. (2015). Furthermore, when operated under strong non-equilibrium conditions where linear response theory breaks down, a theoretical treatment taking into account the full nonlinear properties is needed Haupt et al. (2013); Leijnse et al. (2010); Argüello-Luengo et al. (2015); Sánchez and López (2016). Only recently have these issues been discussed in strongly interacting QD systems Gergs et al. (2015); Seja et al. (2016); Daré and Lombardo (2017).

Another important factor for thermoelectric effects in CCQD systems is the coupling to the leads which is usually treated in the wide-band approximation assuming energy-independent couplings Bruus and Flensberg (2004). However, energy-dependent couplings to the leads occur naturally in many QD systems Waissman et al. (2013); Thierschmann et al. (2015); Bischoff et al. (2015); Keller et al. (2016), and add an important degree of tunability to the system, and is as crucial for the thermoelectric properties Strasberg et al. (2013); Sánchez and Büttiker (2011); Zhang et al. (2015) as it is for Coulomb drag Sánchez et al. (2010); Bischoff et al. (2015); Kaasbjerg and Jauho (2016); Keller et al. (2016); Lim et al. (2016).

In this work, we present a master-equation approach for the calculation of the nonlinear electronic charge and heat currents in interacting QD systems which takes into account the above-mentioned factors. The charge and heat transfer rates produced by electron tunneling processes are obtained with a perturbative $T$-matrix approach Bruus and Flensberg (2004) which allows us to treat sequential and cotunneling processes on equal footing. We resolve the technical challenges associated with the evaluation of the cotunneling rates with an implementation of the often applied regularization scheme Turek and Matveev (2002); Koch et al. (2004) which applies to the general case of energy-dependent lead couplings, applied biases, and temperature gradients in the system.

The main findings and the organization of the paper are as follows. In Sec. II, we introduce the model system of CCQDs. In Sec. III, we present the methodology, and benchmark the approach in Sec. IV by comparing it to the Landauer-Büttiker formalism for transport through a non-interacting single-level QD. In Sec. V, we study nonlinear thermoelectric phenomena in CCQDs. We investigate the energy exchange mediated by the inter-dot Coulomb interaction which among other thermoelectric effects leads to the demon-induced cooling mechanism Strasberg et al. (2013); Koski et al. (2015). Our findings shed light on the limitations imposed by cotunneling processes on the performance of this mechanism. Furthermore, we demonstrate a strongly enhanced performance of the demon-induced cooling effect by tuning the energy-dependence of the lead couplings. In such performance optimization, as we show, cotunneling processes are essential for a quantitative description of the thermoelectric properties. Finally, Sec. VI presents our conclusions, and App. A gives technical details on the cotunneling rates and the regularization procedure.

II Coulomb-coupled QD systems

We consider CCQD systems like the one illustrated in Fig. 1, which can be described by the Hamiltonian

[TABLE]

and consists of a system of CCQDs with Hamiltonian $\hat{H}_{\text{dots}}$ which is coupled to external leads with Hamiltonian $\hat{H}_{\text{leads}}$ by tunnel couplings described by $\hat{H}_{T}$. We denote $\hat{H}_{0}=\hat{H}_{\text{dots}}+\hat{H}_{\text{leads}}$.

For the QD system, we consider a minimal spinless model of inter-dot Coulomb-coupled single-level QDs described by the Hamiltonian

[TABLE]

where $\hat{c}_{\delta}^{\dagger}$ ($\hat{c}_{\delta}^{\phantom{\dagger}}$) creates (annihilates) an electron in QD $\delta$ with energy controlled by gate voltages $\epsilon_{\delta}=-eV_{\delta}$, where $V_{\delta}$ is the gate potential on dot $\delta$, $\hat{n}^{\phantom{\dagger}}_{\delta}=\hat{c}_{\delta}^{\dagger}\hat{c}^{\phantom{\dagger}}_{\delta}$ is the occupation number operator, $U_{\delta\delta^{\prime}}$ is the inter-dot Coulomb interaction, and the summation in the second term is over all QD pairs (specific systems are studied in Secs. IV–V).

The leads are described by non-interacting electron reservoirs, $\hat{H}_{\text{leads}}=\sum_{\ell k}\epsilon_{\ell k}^{\phantom{\dagger}}\hat{c}_{\ell k}^{\dagger}\hat{c}_{\ell k}^{\phantom{\dagger}}$, where $\hat{c}_{\ell k}^{\dagger}$ ($\hat{c}_{\ell k}^{\phantom{\dagger}}$) creates (annihilates) an electron with momentum $k$ and energy $\epsilon_{\ell k}$ in lead $\ell$, which is assumed to be in local equilibrium with temperature $T_{\ell}$ and electrochemical potential $\mu_{\ell}=\mu_{0}-eV_{\ell}$, where $\mu_{0}$ is the equilibrium chemical potential and $V_{\ell}$ is the voltage applied to lead $\ell$. The tunneling Hamiltonian which couples the QD system to the leads is $\hat{H}_{T}=\sum_{\ell k\delta}(t^{\phantom{\dagger}}_{\ell k\delta}\hat{c}^{\dagger}_{\delta}\hat{c}^{\phantom{\dagger}}_{\ell k}+\textit{h.c.}),$ where $t_{\ell k\delta}$ is the tunneling amplitude. We define lead coupling strengths as $\gamma^{\ell}(\epsilon)\equiv 2\pi d_{\ell}(\epsilon)|t_{\ell}(\epsilon)|^{2}$, where $d_{\ell}(\epsilon)$ is the lead density of states. $\gamma^{\ell}(\epsilon)$ is allowed to be energy dependent in contrast to the often applied wide-band approximation.

III Master equation and transport currents

We describe the dynamics and transport in the CCQD system with a Pauli ME where the transitions between the QD states are governed by electron tunneling to and from the leads Timm (2008). The tunneling-induced transition rates are calculated based on a perturbative $T$-matrix approach where the tunneling Hamiltonian is treated as a perturbation to the decoupled QD system and leads. This allows a systematic expansion in the tunnel couplings and the inclusion of high-order processes. However, quantum effects such as tunneling-induced level broadening and level shifts König et al. (1996); Thielmann et al. (2005); Pedersen and Wacker (2005) are not captured by this perturbative approach, which is only valid in the weak coupling regime $\gamma<k_{B}T,U$.

In the absence of tunnel coupling, the states of the decoupled QD system and leads are described by product states of the QD system occupation states $\lvert m\rangle$ with energy $E_{\text{dots},m}=\langle m\rvert\hat{H}_{\text{dots}}\lvert m\rangle$ and the leads $\lvert i\rangle$ with energy $E_{\text{leads},i}=\langle i\rvert\hat{H}_{\text{leads}}\lvert i\rangle$. The non-equilibrium occupations of the QD states are described by probabilities $p_{m}$ (the diagonal components of the reduced density operator of the CCQD system) which are determined by the ME

[TABLE]

where $\Gamma_{mn}$ denotes the tunneling-induced transition rate from QD state $\lvert m\rangle$ to $\lvert n\rangle$. The ME is solved for the steady-state probabilities, $\dot{p}_{m}=0$, in the following. The QD states are given explicitly in Sec. IV and Sec. V for the considered systems.

III.1 Transition rates

The rates for transitions between the QD states are obtained from the generalized Fermi’s golden rule Bruus and Flensberg (2004); Timm (2011)

[TABLE]

where $\Delta_{mn}\equiv E_{\text{dots},n}-E_{\text{dots},m}$, $\rho_{i}$ is the thermal probability of finding the leads in the initial state, the sum is over initial and final states of the leads, and the $T$ matrix obeys

[TABLE]

with $E_{\text{initial}}=E_{\text{dots},m}+E_{\text{leads},i}$, and $\eta$ is a positive infinitesimal.

The lowest-order contribution to the tunneling rates describes single-electron tunneling, or sequential tunneling, processes between the QD system and the leads:

[TABLE]

where Eq. (6) (Eq. (7)) is the sequential rate of tunneling out of, $\rightarrow$, (into, $\leftarrow$) lead $\ell$, thereby changing the state of the QD system from $m$ to $n$, ${f^{\ell}(\epsilon)=[\exp{(\beta_{\ell}(\epsilon-\mu_{\ell}))}+1]^{-1}}$ is the Fermi-Dirac distribution in lead $\ell$, $\bar{f}^{\ell}(\epsilon)=1-f^{\ell}(\epsilon)$, and $\beta_{\ell}=1/(k_{B}T_{\ell})$. The leads are assumed to equilibrate to the Fermi-Dirac distribution in between the tunneling events.

The next-to-leading order terms in the $T$ matrix describe cotunneling processes. In conventional local elastic and inelastic cotunneling processes, a net electron is transferred between two leads attached to the same QD (e.g., System 1 in Fig. 1). Here we also consider (i) nonlocal cotunneling processes Amasha et al. (2013); Kaasbjerg and Jauho (2016) in which a net electron is transferred between leads attached to different QDs, as well as (ii) pair-cotunneling processes where two electrons tunnel into/out of the CCQD system in one coherent process Koch et al. (2006a); Leijnse et al. (2009).

For the thermoelectric effects in focus here, the process of nonlocal cotunneling is important. The (unregularized) rate for nonlocal cotunneling which net transfers an electron out of lead $\ell$ and into lead $\ell^{\prime}$ is given by

[TABLE]

where $v$ ($v^{\prime}$) refers to the virtually occupied intermediate state in the process where an electron initially tunnels from lead $\ell$ and into the QD system (from the QD system and into lead $\ell^{\prime}$). We refer to App. A for the expressions for the remaining cotunneling processes relevant for this study.

A well-known artifact of the cotunneling rates obtained with the $T$-matrix approach is that they formally diverge in the limit $\eta\to 0$. To deal with this divergence different regularization schemes have been proposedTurek and Matveev (2002); Koch et al. (2004); Koller et al. (2010); Timm (2011). Deep inside the Coulomb blockade, the discripancy between the different regularization schemes vanishesKoller et al. (2010). In this work, we apply the by now standard regularization scheme in Ref. Turek and Matveev, 2002, but for future work, a detailed comparison of the charge and heat currents obtained from different regularization schemes could be useful. We denote the regularized rates which enter into Eq. (3) without a tilde. To be explicit, we consider the processes $\Gamma^{\phantom{\ell\rightarrow}}_{mn}\equiv\sum_{\ell}(\Gamma_{mn}^{\ell\leftarrow}+\Gamma_{mn}^{\ell\rightarrow}),\ \Gamma^{\ell\leftarrow}_{mn}\equiv\Gamma^{\small{\overleftarrow{\ell\phantom{{}^{\prime}}}}}_{mn}+\sum_{\ell^{\prime}}(\Gamma^{\small{\overleftarrow{\ell\phantom{{}^{\prime}}}\overrightarrow{\ell^{\prime}}}}_{mn}+\Gamma^{\small{\overleftarrow{\ell\phantom{{}^{\prime}}}\overleftarrow{\ell^{\prime}}}}_{mn}),\ \Gamma^{\ell\rightarrow}_{mn}\equiv\Gamma^{\small{\overrightarrow{\ell\phantom{{}^{\prime}}}}}_{mn}+\sum_{\ell^{\prime}}(\Gamma^{\small{\overrightarrow{\ell\phantom{{}^{\prime}}}\overleftarrow{\ell^{\prime}}}}_{mn}+\Gamma^{\small{\overrightarrow{\ell\phantom{{}^{\prime}}}\overrightarrow{\ell^{\prime}}}}_{mn}).$ A numerical procedure for the regularization is outlined in App. A.

III.2 Charge and heat currents

The steady-state transport currents can be obtained from the occupation probabilities. The electric current going into lead $\ell$ is

[TABLE]

where $\hat{n}_{\ell k}^{\phantom{\dagger}}=\hat{c}_{\ell k}^{\dagger}\hat{c}_{\ell k}^{\phantom{\dagger}}$, $p_{m}$ is calculated from Eq. (3), and the rightmost form expresses the electric current in terms of the total rate of electrons tunneling into lead $\ell$, minus the total rate of electrons tunneling out of lead $\ell$FN (2).

The heat current going into lead $\ell$ is Leijnse et al. (2010); Ludovico et al. (2016); Benenti et al. (2016)

[TABLE]

where the rightmost form expresses the heat current in terms of heat rates $W$ (using a similar notation as for the tunneling rates).

The sequential-tunneling heat rate in lead $\ell$ is calculated as the tunneling rate multiplied by the energy of the tunneling electron relative to the chemical potential in the lead,

[TABLE]

where the indices follow the notation of the tunneling rates, however, the additional first subscript $\ell$ refers to the lead in which the heat rate is calculated.

Analogously, the cotunneling heat rates into/out of the leads are calculated a posteriori by multiplying the integrand in the cotunneling rate by the energy of the tunneling electron relative to the chemical potential of the lead. For example, for the nonlocal cotunneling process between lead $\ell$ and $\ell^{\prime}$, the heat rate in lead $\ell$ reads

[TABLE]

with the heat rate in lead $\ell^{\prime}$, $\tilde{W}^{\small{\overrightarrow{\ell\phantom{{}^{\prime}}}\overleftarrow{\ell^{\prime}}}}_{\ell^{\prime},mn}$, given as above but with $(\epsilon-\mu_{\ell})$ replaced by $(\epsilon-\Delta_{mn}-\mu_{\ell^{\prime}})$. The remaining cotunneling heat rates follow similarly.

Whereas the calculation of charge currents involves the electron-tunneling rates which enter the ME (3), and therefore does not require any additional steps once the ME has been set up and solved, the heat currents must be calculated via the heat tunneling rates in a post-processing step, similar to more rigorous density-matrix treatments. Gergs et al. (2015)

IV Comparison to the Landauer-Büttiker formalism

In this section, we benchmark the approach by comparing the charge and heat currents in a spinless non-interacting single-level QD system with those obtained from the Landauer-Büttiker (LB) formalism (see Ref. Koch et al., 2006b for a comparison of the electric current in the case of equal temperatures in the leads). For non-interacting systems the LB result is exact. However, for the thermoelectric effects discussed in Sec. V which require the presence of strong Coulomb interaction, an alternative method such as the ME approach is needed.

We consider a single-level QD coupled to two leads $\ell\in\{A,B\}$ (such as System 1 in Fig. 1 when tunnel- and Coulomb-decoupled from System 2). For simplicity, we assume energy-independent lead couplings $\gamma^{\ell}(\epsilon)=\gamma^{\ell}$ in this case. The Hamiltonian of the QD reduces to

[TABLE]

with states labeled by the occupancy, $\lvert n_{1}\rangle\in\{\lvert 0\rangle,\lvert 1\rangle\}$.

In the LB formalism, the electric current and heat current going into lead $A$ are given by Haug and Jauho (2008); Argüello-Luengo et al. (2015),

[TABLE]

and

[TABLE]

respectively. For a non-interacting single-level QD the transmission function $T(\epsilon)$ is

[TABLE]

where $\gamma=\gamma^{A}+\gamma^{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 approaches with a finite bias and temperature difference ($T_{B}=2T_{A}\equiv 2T$) between the leads are plotted in Figs. 2(a) and 2(b) as a function of the gate voltage for two different lead coupling strengths. To demonstrate the importance of cotunneling processes, we have included ME results based on sequential tunneling only (black dotted curves) which do not depend on $\gamma^{\ell}$ in the units shown, as well as sequential plus cotunneling (dashed curves). The results based purely on sequential tunneling differ significantly from the LB results unless $\gamma^{\ell}\ll k_{B}T$. However, for $\gamma^{\ell}<k_{B}T$, the ME results with cotunneling are in excellent agreement with the LB formalism. For $\gamma^{\ell}>k_{B}T$ which is outside the regime of validity of the ME approach, the two approaches deviate, as expected.

In the following discussion of thermoelectric effects, the heat current is of particular interest. As seen in Fig. 2(b), when the dot level is above the electrochemical potential in lead $A$, the heat current becomes negative (for sufficiently small lead coupling strength). In this case, electrons above the electrochemical potential tunnel out of the lead and thereby cool the lead [cf. Eq. (10)]. Such cooling mechanisms due to energy-selective tunneling have been confirmed experimentally in metallic QD systems Feshchenko et al. (2014); Koski et al. (2015). The energy-selective tunneling gives rise to an asymmetry in the energy dissipation between the source and drain leads which was recently observed in molecular junctions Lee et al. (2013).

V Thermoelectric effects in Coulomb-coupled QDs

In the remaining part of the paper, we study the thermoelectric properties of the system illustrated in Fig. 1, i.e. two single-level QDs with QD1 tunnel-coupled to leads $A$ and $B$ and QD2 tunnel-coupled to lead $C$. The CCQD system is described by the Hamiltonian

[TABLE]

where we have used the simplified notation $U_{12}\equiv U$, and the occupation states are $\lvert m\rangle=\lvert n_{1}n_{2}\rangle\in\{\lvert 00\rangle,\lvert 10\rangle,\lvert 01\rangle,\lvert 11\rangle\}$. We consider situations where a source-drain bias $V$ is applied to System 1, $\mu_{A}=\mu_{0}+eV/2$, $\mu_{B}=\mu_{0}-eV/2$ (we set $\mu_{0}=0$ as reference).

As pointed out above, we here allow for energy-dependent lead couplings. For bias voltages and temperature differences small compared to the energy scale at which the lead couplings vary, it suffices to consider the expansion of the lead couplings around their value at $\mu_{0}$FN (1),

[TABLE]

where $\gamma_{0}^{\ell}=\gamma^{\ell}(\mu_{0})$, $\partial\gamma^{\ell}\equiv\tfrac{\partial\gamma^{\ell}(\epsilon)}{\partial\epsilon}|_{\epsilon=\mu_{0}}$.

V.1 Current and energy exchange

In Fig. 3(a) we show the electric current through QD1, $I\equiv I_{A}=-I_{B}$, at low temperature $k_{B}T_{\ell}=10^{-2}U$ (for illustrative convenience) and bias $eV=0.3U$ as a function of gate detuning $V_{2}-V_{1}$ and total gating $V_{1}+V_{2}$ in the vicinity of the honeycomb vertex of the stability diagram van der Wiel et al. (2002). Here, we initially assume energy-independent lead couplings which is sufficient to get an overall understanding of the behavior of the system. The large current near the degeneracy lines defined by $\Delta_{00,10}=0$ and $\Delta_{01,11}=0$ is due to sequential tunneling processes. Away from these degeneracy lines where sequential tunneling is exponentially suppressed, cotunneling processes give rise to a weak background current. At the degeneracy line $\Delta_{10,01}=0$ connecting the two triple points at $V_{1}=V_{2}=0,U$, respectively, nonlocal cotunneling processes are responsible for the enhanced cotunneling current.

The heat currents which accompany the electric current are shown in Figs. 3(b)–3(d) for different temperatures in the leads. Figure 3(b) shows the heat current in lead $A$ for $k_{B}T_{\ell}=0.1U$. Along the degeneracy lines where $\Delta_{00,10}=0$ and $\Delta_{01,11}=0$ and only the occupation of QD1 fluctuates, the heat current shows a behavior similar to the one in Fig. 2(b) for a single-level QD. However, at the center of the stability diagram, Coulomb-mediated energy exchange due to the strong Coulomb interaction between the QDs becomes significant. This manifests itself in a cooling of System 1 inside the region bounded by the solid lines at the center of Fig. 3(b) (notice that the color scale is dominated by the heat current with larger magnitude outside this region). From the heat current in lead $C$ shown in Fig. 3(c), the cooling of System 1 is seen to be at the cost of heating System 2. This Coulomb-mediated energy exchange between the two QD systems occurs in spite of the fact that no electrons are exchanged, and is the driving force behind demon-induced cooling Strasberg et al. (2013); Koski et al. (2015), energy harvesting Sánchez and Büttiker (2011); Sothmann et al. (2014); Thierschmann et al. (2015, 2016), and Coulomb drag Sánchez et al. (2010); Kaasbjerg and Jauho (2016).

A simple analytical result for the energy exchange can be found when considering sequential tunneling processes only (indicated by the superscript $s$). In this case, the total heat currents in System 1, $J_{1}^{s}\equiv J_{A}^{s}+J_{B}^{s}$, and System 2, $J_{2}^{s}\equiv J_{C}^{s}$, become Sánchez and Büttiker (2011)

[TABLE]

where $\Gamma^{-}\equiv\Gamma_{00,01}\Gamma_{01,11}\Gamma_{11,10}\Gamma_{10,00}$, $\Gamma^{+}\equiv\Gamma_{00,10}\Gamma_{10,11}\Gamma_{11,01}\Gamma_{01,00}$. The factor $\tau^{s}$ depends on the various sequential tunneling rates, however, is merely a normalization factor and is not reproduced here. The first two terms proportional to $U$ in Eq. (19) describe the energy exchange, whereas the last term in Eq. (19a) describes the contribution from Joule heating in System 1. The direction of the energy transfer is determined by the sign of $\Gamma^{-}-\Gamma^{+}$. It is therefore convenient to consider the ratio

[TABLE]

which describes whether energy is transferred from System 1 to 2 ($\Gamma^{-}/\Gamma^{+}>1$) or vice versa ($\Gamma^{-}/\Gamma^{+}<1$) Sánchez and Büttiker (2012). On the right-hand side of (20), we have taken $\beta_{A/B}=\beta_{1}$ and $\beta_{C}=\beta_{2}$, and expressed the ratio in terms of an exponential factor, which depends on the temperature in System 1 and System 2, and

[TABLE]

which depends on the temperature and bias in System 1 only. The subscript 0 (1) in Eq. (21) indicates that the corresponding function is evaluated at $\Delta_{00,10}$ ($\Delta_{01,11}$) [see Eqs. (6)–(7)].

The exponential factor in (20) shows that a temperature gradient between the two QD systems can generate a net heat flow from the hot to the cold system. This is the mechanism behind the heat engine studied in Ref. Sánchez and Büttiker, 2011. On the other hand, a closer inspection of the $\Omega$ factor reveals that it is, in fact, possible to generate a net heat flow in the opposite direction, i.e. from the cold to the hot system, and this is the cause of the negative heat current at the center of Fig. 3(b). This so-called demon-induced cooling effect will be discussed further in Sec. V.2 below.

When the applied bias and temperature are small compared to the inter-dot Coulomb interaction, $eV,k_{B}T\ll U$, cotunneling processes start to dominate the heat currents. This is demonstrated in Fig. 3(d) which shows the heat current $J_{C}$ as a function of temperature for the two different gate tunings marked with symbols in Fig. 3(c). Considering sequential tunneling only (dashed curves), the heat current is quenched at $k_{B}T\ll U$ as $\Gamma_{01,11}$ and $\Gamma_{10,00}$ in $\Gamma^{-}$ become exponentially suppressed. This can also be understood from the illustration in Fig. 4(a) which shows the sequence of sequential tunneling processes corresponding to $\Gamma^{-}$. However, nonlocal cotunneling processes allow the system to fluctuate between the two states $10\leftrightarrow 01$, as illustrated in Fig. 4(b), and thereby transfer heat between the systems. The nonlocal cotunneling channel is open for $|\Delta_{01,10}|\lesssim\text{max}\{|eV/2|,k_{B}T\}$, and the associated heat current is thus also suppressed at low temperature when $\Delta_{01,10}\neq 0$ as illustrated by the blue curve (triangle) in Fig. 3(d). For zero detuning $\Delta_{01,10}=0$ (circle), the nonlocal cotunneling rates, and hence also the heat current, saturate at $k_{B}T\ll eV$. In Sec. V.2, we demonstrate that nonlocal cotunneling processes have a significant effect on the demon-induced cooling mechanism.

V.2 Demon-induced cooling

The effect of cooling System 1 at the cost of heating System 2 has recently been discussed in context of a Maxwell’s demon where System 2 plays the role of the demon which performs the necessary feedback to cool System 1 Strasberg et al. (2013); Koski et al. (2015). To maximize the achievable cooling power for refrigeration purposesKutvonen et al. (2016), large tunneling rates, $\gamma^{\ell}(\epsilon)\sim k_{B}T,U$, are essential [cf. Eq. (19)]. However, large tunneling rates increase the contribution from higher-order tunneling processes, thus emphasizing the importance of including cotunneling processes in quantitative analyses of the cooling power.

In the following, we consider the case of uniform temperature $T_{\ell}\equiv T$ whereby the exponential factor in (20) becomes unity. This allows us to focus on the $\Omega$ factor in the optimization of the performance. Equation (19) shows that the cooling mechanism is governed by $\Gamma^{-}$ since, as illustrated in Fig. 4(a), in a full sequential cycle an amount of energy $U$ is transferred from System 1 to System 2 thereby cooling System 1. In the following, we discuss how to increase the cooling power by maximizing the success rate for completing the cooling cycle in Fig. 4(a).

V.2.1 Cotunneling limitations

Although the cycle of nonlocal cotunneling processes illustrated in Fig. 4(b) gives the same net transfer of electrons as the sequential tunneling cycle in Fig. 4(a), the net energy transfer is different for the two cases. As illustrated, in a cotunneling process also electrons below (above) the electrochemical potential can tunnel out of lead $A$ (into lead $B$), and thus reduce the demon-induced cooling effect.

In Fig. 4(c), we show the heat current $J_{1}$ together with its individual contributions from sequential ($J_{1}^{s}$) and cotunneling ($J_{1}^{c}$) processes. Overall, System 1 cools at low bias, while at higher bias, Joule heating becomes dominant. The minimum in $J_{1}$ as a function of bias voltage is referred to as the maximum cooling power, $J_{1,\max}\equiv\min J_{1}(V)$. As the figure shows, cotunneling reduces the maximum cooling power.

Figure 5 shows how the maximum cooling power $J_{1,\max}$ scales with the lead coupling strengths. As the figure demonstrates, the rates must satisfy $\gamma^{C}>\gamma^{A/B}$ to ensure that System 2 acts sufficiently fast to perform the desired feedback such that the cooling cycle in Fig. 4(a) is completed when an electron tunnels between lead $A$ and $B$ Kutvonen et al. (2016). In the region of large cooling power, cotunneling processes start to become important, and hence there is a trade-off between sequential tunneling which improves the cooling effect, and nonlocal cotunneling which limits the effect. In addition, the area in the lead coupling parameter space where refrigeration is possible is also reduced when cotunneling is included.

V.2.2 Performance boosting

Here we demonstrate that energy-dependent lead couplings can enhance the demon-induced cooling power significantly. We restrict the discussion to lead couplings with a linear energy dependence [cf. Eq. (18)].

By inspecting the $\Omega$ factor in Eq. (21), we find that for $\mu_{A}>\mu_{B}$, the configuration illustrated in the inset of Fig. 6 where $\gamma^{A}_{0}$, $\gamma^{B}_{1}$ are reduced compared to $\gamma^{A}_{1}$, $\gamma^{B}_{0}$, boosts the $\Omega$ factor (and thereby $\Gamma^{-}/\Gamma^{+}$). This results in an enhancement of the cooling power by suppressing direct tunneling between lead $A$ and $B$ via two sequential tunneling processes which contributes to Joule heating, while at the same time promoting the processes of the cooling cycle in Fig. 4(a).

In Fig. 6 we show the maximum cooling power as a function of temperature for different situations for the energy dependence of the lead couplings, from the top (black) curve showing the result for energy independent lead couplings, to increasing energy dependence, i.e. increasing $\lvert\partial\gamma^{A/B}\rvert$, towards the bottom (light blue) curve. When tuning the energy dependence of the lead couplings, a significant enhancement of the cooling power is achieved. Again, the effect of cotunneling processes is to reduce the attainable cooling power (solid lines) relative to the cooling power obtained when only considering sequential tunneling processes (dashed lines).

VI Conclusions

In summary, we have studied thermoelectric effects in CCQD systems with a $T$-matrix based master-equation approach for the calculation of charge and heat currents. Importantly, our method (i) treats incoherent sequential tunneling processes and coherent cotunneling processes on equal footing, and (ii) can account for energy-dependent tunnel couplings to the leads. Both are essential for quantitative predictions and optimization of the thermoelectric properties of CCQDs.

To benchmark the master-equation method, we considered a non-interacting single-level QD coupled to source and drain leads for which the Landauer-Büttiker formalism is exact. In the regime of validity of our method, i.e. small tunnel couplings to the leads, $\gamma<k_{B}T$, we demonstrated excellent agreement with the results from the Landauer-Büttiker method when cotunneling processes are included in the master equation.

Furthermore, we studied the effect of cotunneling processes and energy-dependent lead couplings on the thermoelectric properties of a CCQD system consisting of two QDs exhibiting a Maxwell’s demon-like cooling mechanism Strasberg et al. (2013); Koski et al. (2015). First of all, we showed that cotunneling processes reduce the cooling effect since cotunneling processes do not share the delicate energy selectivity inherent to sequential tunneling processes. This results in a significant reduction of the achievable cooling power compared to the sequential tunneling result when the lead couplings are increased to maximize the cooling power from sequential tunneling processes. Secondly, we demonstrated that it is possible to boost the cooling power significantly via other means by introducing energy-dependent lead couplings and properly tuning their energy dependence. In this case, we showed that cotunneling still reduces the cooling power significantly, thus emphasizing the importance of cotunneling processes in quantitative analyses.

Applying the methodology to other mesoscopic systems allows for testing of new thermoelectric device ideas beyond sequential tunneling estimates, as well as for improved comparison with experiments.

Acknowledgements.

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’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 713683 (COFUNDfellowsDTU). The Center for Nanostructured Graphene (CNG) is sponsored by the Danish Research Foundation, Project DNRF103.

Appendix A Cotunneling rates and regularization procedure

The rate for elastic cotunneling through a single-level QD is given by

[TABLE]

where $v$ refers to the virtually occupied intermediate state created in the process where an initially empty level is filled ($+\epsilon$) or an initially filled level is emptied ($-\epsilon$).

In pair-cotunneling processes, two electrons tunnel simultaneously out of (into) the QD system and into (out of) the leads $\ell$ and $\ell^{\prime}$. The rate for such processes takes the form

[TABLE]

where $v$ ($v^{\prime}$) refers to the virtually occupied intermediate state in a process where an electron initially tunnels from the QD system and into lead $\ell$ ($\ell^{\prime}$). Similarly,

[TABLE]

where $v$ ($v^{\prime}$) refer to the virtually occupied intermediate state in a process where an electron initially tunnels from lead $\ell^{\prime}$ ($\ell$) and into the QD system.

The bare cotunneling rates are formally divergent in the limit $\eta\to 0$. The divergence stems from factors involving $|x+i\eta|^{-2}$, $x,\eta\in\mathbb{R}$. Using that Turek and Matveev (2002)

[TABLE]

where $\mathcal{P}$ denotes the principle value, we can identify the divergent contributions, e.g. from Eq. (III.1)

[TABLE]

where $\Gamma_{mn}^{\small{\overrightarrow{\ell\phantom{{}^{\prime}}}\overleftarrow{\ell^{\prime}}}}$ denotes the regularized cotunneling rate, and we have used that the cross-terms from the absolute squared in Eq. (III.1) do not contribute to any divergences. The divergent contribution is proportional to products of two sequential tunneling rates. These correspond to two energy-conserving (sequential) transitions which can be identified with the intermediate processes in the cotunneling process. The sum is over the possible sequences of intermediate transitions. Similarly, for the cotunneling heat rates, e.g. Eq. (III.2)

[TABLE]

or the corresponding heat rate in lead $\ell^{\prime}$

[TABLE]

We apply the regularization scheme in Ref. Turek and Matveev, 2002 and subtract these terms scaling as $\eta^{-1}$.

In the case of identical temperatures in the leads, using the identity $f(\epsilon_{1})[1-f(\epsilon_{2})]=n(\epsilon_{1}-\epsilon_{2})[f(\epsilon_{2})-f(\epsilon_{1})]$, where $f(\epsilon)$ is the Fermi-Dirac distribution and $n(\epsilon)$ is the Bose-Einstein distribution, the cotunneling rates can be written in the form

[TABLE]

where $P(\epsilon)$ is assumed to be a polynomial, $P(\epsilon)=\sum_{i=0}^{n}c_{i}\epsilon^{n}$, of maximum order $n=2$ for $k_{1}-k_{2}\neq 0$ and $n=4$ for $k_{1}-k_{2}=0$ to ensure that the result below is well-defined. The derivation is in line with the one in Ref. Kaasbjerg and Jauho, 2016, and the integral becomes

[TABLE]

where

[TABLE]

with $\psi$ ($\psi_{1}$) being the digamma (trigamma) function, and

[TABLE]

The term $\mathcal{O}(\eta^{-1})$ is omitted by regularization before letting $\eta\to 0$. For $k_{B}T<\gamma$ (outside the regime of validity), the failure of the approach is seen as a logarithmic divergence of the digamma functions near the degeneracy points.

In studies of thermoelectric effects where different lead temperatures as well as more general energy-dependence of the lead couplings become relevant, one must turn to a numerical procedure. In this case, we evaluate the cotunneling integrals numerically with a small but finite $\eta$, and subsequently subtract contributions of order $\eta^{-1}$ as shown in e.g. Eqs. (26)–(28). In particular, we have applied the numerical procedure in Fig. 2 and Fig. 6, and stated the values of $\eta$ in the figure caption.