Thermoelectric properties of an interacting quantum dot-based heat   engine

Paolo Andrea Erdman; Francesco Mazza; Riccardo Bosisio; Giuliano; Benenti; Rosario Fazio; Fabio Taddei

arXiv:1702.06042·cond-mat.mes-hall·June 29, 2017

Thermoelectric properties of an interacting quantum dot-based heat engine

Paolo Andrea Erdman, Francesco Mazza, Riccardo Bosisio, Giuliano, Benenti, Rosario Fazio, Fabio Taddei

PDF

TL;DR

This paper analyzes the thermoelectric performance of an interacting quantum dot heat engine, deriving analytical expressions for transport coefficients and optimizing conditions for maximum power and efficiency.

Contribution

It provides new analytical formulas for transport properties of multi-level quantum dots in the quantum limit, considering interactions and level degeneracy.

Findings

01

Optimal conditions for maximum power and efficiency identified.

02

Energy level degeneracy enhances output power.

03

Breakdown of Onsager reciprocity under large biases observed.

Abstract

We study the thermoelectric properties and heat-to-work conversion performance of an interacting, multi-level quantum dot (QD) weakly coupled to electronic reservoirs. We focus on the sequential tunneling regime. The dynamics of the charge in the QD is studied by means of master equations for the probabilities of occupation. From here we compute the charge and heat currents in the linear response regime. Assuming a generic multi-terminal setup, and for low temperatures (quantum limit), we obtain analytical expressions for the transport coefficients which account for the interplay between interactions (charging energy) and level quantization. In the case of systems with two and three terminals we derive formulas for the power factor Q and the figure of merit ZT for a QD-based heat engine, identifying optimal working conditions which maximize output power and efficiency of heat-to-work…

Figures16

Click any figure to enlarge with its caption.

Equations274

f_{α} (E) = [1 + exp (\frac{E - μ _{α}}{k _{B} T _{α}})]^{- 1},

f_{α} (E) = [1 + exp (\frac{E - μ _{α}}{k _{B} T _{α}})]^{- 1},

E_{p} + U (N) = E^{fin} (N) + U (N - 1) .

E_{p} + U (N) = E^{fin} (N) + U (N - 1) .

E^{in} (N) + U (N) = E_{p} + U (N + 1) .

E^{in} (N) + U (N) = E_{p} + U (N + 1) .

\frac{\partial}{\partial t} P ({n_{i}}) = - p α \sum δ_{n_{p}, 0} P ({n_{i}}) Γ_{α} (p) f_{α} (E^{in} (N)) - p α \sum δ_{n_{p}, 1} P ({n_{i}}) Γ_{α} (p) [1 - f_{α} (E^{fin} (N))] + p α \sum δ_{n_{p}, 0} P ({n_{i}}, n_{p} = 1) Γ_{α} (p) [1 - f_{α} (E^{fin} (N + 1))] + p α \sum δ_{n_{p}, 1} P ({n_{i}}, n_{p} = 0) Γ_{α} (p) f_{α} (E^{in} (N - 1)),

\frac{\partial}{\partial t} P ({n_{i}}) = - p α \sum δ_{n_{p}, 0} P ({n_{i}}) Γ_{α} (p) f_{α} (E^{in} (N)) - p α \sum δ_{n_{p}, 1} P ({n_{i}}) Γ_{α} (p) [1 - f_{α} (E^{fin} (N))] + p α \sum δ_{n_{p}, 0} P ({n_{i}}, n_{p} = 1) Γ_{α} (p) [1 - f_{α} (E^{fin} (N + 1))] + p α \sum δ_{n_{p}, 1} P ({n_{i}}, n_{p} = 0) Γ_{α} (p) f_{α} (E^{in} (N - 1)),

P ({n_{i}}, n_{p} = 1) = P ({n_{1}, \dots, n_{p - 1}, 1, n_{p + 1}, \dots})

P ({n_{i}}, n_{p} = 1) = P ({n_{1}, \dots, n_{p - 1}, 1, n_{p + 1}, \dots})

P ({n_{i}}, n_{p} = 0) = P ({n_{1}, \dots, n_{p - 1}, 0, n_{p + 1}, \dots})

P ({n_{i}}, n_{p} = 0) = P ({n_{1}, \dots, n_{p - 1}, 0, n_{p + 1}, \dots})

\sum_{{n_{i}}} P ({n_{i}}) = 1

\sum_{{n_{i}}} P ({n_{i}}) = 1

J^{c}_{\alpha}=e\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))\\ -\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]\Big{\}},

J^{c}_{\alpha}=e\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))\\ -\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]\Big{\}},

J^{u}_{\alpha}=\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))E^{\text{in}}(N)-\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]E^{\text{fin}}(N)\Big{\}},

J^{u}_{\alpha}=\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))E^{\text{in}}(N)-\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]E^{\text{fin}}(N)\Big{\}},

J^{h}_{\alpha}=\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))\left[E^{\text{in}}(N)-\mu_{\alpha}\right]-\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]\left[E^{\text{fin}}(N)-\mu_{\alpha}\right]\Big{\}}.

J^{h}_{\alpha}=\sum\limits_{p=1}^{\infty}\sum\limits_{\{n_{i}\}}P(\{n_{i}\})\Gamma_{\alpha}(p)\Big{\{}\delta_{n_{p},0}f_{\alpha}(E^{\text{in}}(N))\left[E^{\text{in}}(N)-\mu_{\alpha}\right]-\delta_{n_{p},1}[1-f_{\alpha}(E^{\text{fin}}(N))]\left[E^{\text{fin}}(N)-\mu_{\alpha}\right]\Big{\}}.

p \sum (δ_{n_{p}, 1} - δ_{n_{p}, 0}) [P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p}] = 0,

p \sum (δ_{n_{p}, 1} - δ_{n_{p}, 0}) [P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p}] = 0,

A_{N, p} = α \sum Γ_{α} (p) f_{α} (E^{in} (N))

A_{N, p} = α \sum Γ_{α} (p) f_{α} (E^{in} (N))

B_{N, p} = α \sum Γ_{α} (p) [1 - f_{α} (E^{in} (N))] = α \sum Γ_{α} (p) - A_{N, p} .

B_{N, p} = α \sum Γ_{α} (p) [1 - f_{α} (E^{in} (N))] = α \sum Γ_{α} (p) - A_{N, p} .

E^{fin} (N + 1) = E^{in} (N),

E^{fin} (N + 1) = E^{in} (N),

P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p} = 0

P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p} = 0

α, {n_{i}}_{i \neq = p} \sum P ({n_{i}}, n_{p} = 0) Γ_{α} (p) f_{α} (E^{in} (N)) .

α, {n_{i}}_{i \neq = p} \sum P ({n_{i}}, n_{p} = 0) Γ_{α} (p) f_{α} (E^{in} (N)) .

α, {n_{i}}_{i \neq = p} \sum P ({n_{i}}, n_{p} = 1) Γ_{α} (p) [1 - f_{α} (E^{fin} (N))] .

α, {n_{i}}_{i \neq = p} \sum P ({n_{i}}, n_{p} = 1) Γ_{α} (p) [1 - f_{α} (E^{fin} (N))] .

{n_{i}}_{i \neq = p} \sum [P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p}] = 0.

{n_{i}}_{i \neq = p} \sum [P ({n_{i}}, n_{p} = 0) A_{N, p} - P ({n_{i}}, n_{p} = 1) B_{N, p}] = 0.

P = α = 1 \sum N J_{α}^{h} .

P = α = 1 \sum N J_{α}^{h} .

η = \frac{P}{α ^{'} \sum J _{α^{'}}^{h}},

η = \frac{P}{α ^{'} \sum J _{α^{'}}^{h}},

P_{max} = \frac{1}{4} Q Δ T^{2},

P_{max} = \frac{1}{4} Q Δ T^{2},

η (P_{max}) = \frac{η _{C}}{2} \frac{Z T}{Z T + 2},

η_{max} = η_{C} \frac{1 + Z T - 1}{1 + Z T + 1},

P ({n_{i}}) = P_{eq} ({n_{i}}) [1 + ψ ({n_{i}})],

P ({n_{i}}) = P_{eq} ({n_{i}}) [1 + ψ ({n_{i}})],

P_{eq} ({n_{i}}) = \frac{1}{Z} exp [- \frac{1}{k _{B} T} (p = 1 \sum \infty E_{p} n_{p} + U (N) - μ N)]

P_{eq} ({n_{i}}) = \frac{1}{Z} exp [- \frac{1}{k _{B} T} (p = 1 \sum \infty E_{p} n_{p} + U (N) - μ N)]

Z = {n_{i}} \sum exp [- \frac{1}{k _{B} T} (p = 1 \sum \infty E_{p} n_{p} + U (N) - μ N)] .

Z = {n_{i}} \sum exp [- \frac{1}{k _{B} T} (p = 1 \sum \infty E_{p} n_{p} + U (N) - μ N)] .

{n_{i}}_{i \neq = p} \sum P_{eq} ({n_{i}}, n_{p} = 0) f (E_{p} + (2 N + 1) E_{C}) \times α \sum Γ_{α} (p) {ψ ({n_{i}}, n_{p} = 0) - ψ ({n_{i}}, n_{p} = 1) + \frac{1}{k _{B} T} [(E_{p} + (2 N + 1) E_{C} - μ) \frac{Δ T _{α}}{T} + Δ μ_{α}]} = 0,

{n_{i}}_{i \neq = p} \sum P_{eq} ({n_{i}}, n_{p} = 0) f (E_{p} + (2 N + 1) E_{C}) \times α \sum Γ_{α} (p) {ψ ({n_{i}}, n_{p} = 0) - ψ ({n_{i}}, n_{p} = 1) + \frac{1}{k _{B} T} [(E_{p} + (2 N + 1) E_{C} - μ) \frac{Δ T _{α}}{T} + Δ μ_{α}]} = 0,

J_{α}^{c} = \frac{e}{k _{B} T} p = 1 \sum \infty N = 1 \sum \infty P_{eq} (N) F_{eq} (E_{p} ∣ N) [1 - f (ϵ (N, p))] β \sum \frac{Γ _{α} ( p ) Γ _{β} ( p )}{Γ _{tot} ( p )} [(ϵ (N, p) - μ) \frac{Δ T _{α} - Δ T _{β}}{T} + (Δ μ_{α} - Δ μ_{β})],

J_{α}^{c} = \frac{e}{k _{B} T} p = 1 \sum \infty N = 1 \sum \infty P_{eq} (N) F_{eq} (E_{p} ∣ N) [1 - f (ϵ (N, p))] β \sum \frac{Γ _{α} ( p ) Γ _{β} ( p )}{Γ _{tot} ( p )} [(ϵ (N, p) - μ) \frac{Δ T _{α} - Δ T _{β}}{T} + (Δ μ_{α} - Δ μ_{β})],

ϵ (N, p) = E_{p} + U (N) - U (N - 1) = E_{p} + (2 N - 1) E_{C}

ϵ (N, p) = E_{p} + U (N) - U (N - 1) = E_{p} + (2 N - 1) E_{C}

P_{eq} (N) \equiv {n_{i}} \sum P_{eq} ({n_{i}}) δ_{\sum n_{i}, N}

P_{eq} (N) \equiv {n_{i}} \sum P_{eq} ({n_{i}}) δ_{\sum n_{i}, N}

F (E_{p} ∣ N) \equiv \frac{P _{eq} ( E _{p} \cap N )}{P _{eq} ( N )} = \frac{{ n _{i} } \sum P _{eq} ({ n _{i} }) δ _{n_{p}, 1} δ _{\sum n_{i}, N}}{P _{eq} ( N )}

F (E_{p} ∣ N) \equiv \frac{P _{eq} ( E _{p} \cap N )}{P _{eq} ( N )} = \frac{{ n _{i} } \sum P _{eq} ({ n _{i} }) δ _{n_{p}, 1} δ _{\sum n_{i}, N}}{P _{eq} ( N )}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Thermoelectric properties of an interacting quantum dot-based heat engine

Paolo Andrea Erdman

NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy

[email protected]

Francesco Mazza

NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy

Riccardo Bosisio

NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy

Giuliano Benenti

Center for Nonlinear and Complex Systems, Dipartimento di Scienza e Alta Tecnologia, Università degli Studi dell’Insubria, via Valleggio 11, 22100 Como, Italy

Istituto Nazionale di Fisica Nucleare, Sezione di Milano, via Celoria 16, 20133 Milano, Italy

NEST, Istituto Nanoscienze-CNR, I-56126 Pisa, Italy

Rosario Fazio

ICTP, Strada Costiera 11, I-34151 Trieste, Italy

NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy

Fabio Taddei

NEST, Istituto Nanoscienze-CNR and Scuola Normale Superiore, I-56126 Pisa, Italy

Abstract

We study the thermoelectric properties and heat-to-work conversion performance of an interacting, multi-level quantum dot (QD) weakly coupled to electronic reservoirs. We focus on the sequential tunneling regime. The dynamics of the charge in the QD is studied by means of master equations for the probabilities of occupation. From here we compute the charge and heat currents in the linear response regime. Assuming a generic multi-terminal setup, and for low temperatures (quantum limit), we obtain analytical expressions for the transport coefficients which account for the interplay between interactions (charging energy) and level quantization. In the case of systems with two and three terminals we derive formulas for the power factor $Q$ and the figure of merit $ZT$ for a QD-based heat engine, identifying optimal working conditions which maximize output power and efficiency of heat-to-work conversion. Beyond the linear response we concentrate on the two-terminal setup. We first study the thermoelectric non-linear coefficients assessing the consequences of large temperature and voltage biases, focusing on the breakdown of the Onsager reciprocal relation between thermopower and Peltier coefficient. We then investigate the conditions which optimize the performance of a heat engine, finding that in the quantum limit output power and efficiency at maximum power can almost be simultaneously maximized by choosing appropriate values of electrochemical potential and bias voltage. At last we study how energy level degeneracy can increase the output power.

pacs:

72.20.Pa,73.23.-b

I Introduction

The study of thermoelectric effects in nanostructures bib:dresselhaus ; bib:vineis ; bib:snyder ; PRreview is attracting increasing interest. Heat-to-work conversion based on thermoelectricity promises an enhanced efficiency as a consequence of the reduction of the phonon contribution to thermal conductance in disordered nanostructures Hochbaum and of the “energy filtering” effect that can result from confinement and quantum effects bib:hicks ; bib:mahan . In particular, an increase of the electron contribution to the figure of merit $ZT$ (which controls the maximum efficiency and the efficiency at maximum power) is possible if one can “filter” the electrons participating in the transport to a narrow energy range bib:mahan .

A heat engine composed of a quantum dot (QD) is a paradigmatic example, since it is characterized by a spectrum of discrete levels which maximizes energy filtering. The thermoelectric properties of QD systems bib:zianni_2007 ; bib:zianni_2008 ; bib:beenakker1 ; bib:beenakker2 ; rejec2012 ; muralidharan2013 ; jacquet2009 ; wohlman2012 ; bjorn2013 ; kennes2013 ; billings2010 ; kubala2006 ; koch2004 ; andreev2001 ; dutt2013 ; costi2010 ; turek2002 ; fazio2001 ; jukka2008 ; mani2011 ; sanchezPRB2013 ; Sanchez2014 ; zimbovskaya2016 and the performance of QD-based heat engines weymann2013 ; donsa2014 ; sothmann2012 ; liu2011 ; kuo2013 ; kuo2010 ; wierzbicki2010 ; rossello2017 ; perroni2016 ; szukiewicz2016 ; Mazza2014 ; kuo2016 ; karol2006 ; murphy2008 ; kuo2009 ; dubi2009 ; imry2010 ; karlstrom2011 ; trocha2012 ; tseng2015 ; taylor2015 ; esposito2009 ; swirkowicz2009 ; leijnse2010 ; nakpathomkun2010 ; buttiker2011 ; muralindharan2012 ; buttiker2013 ; de2016 ; agarwal2014 ; zianni2010 ; liu2010 has been studied theoretically by a number of authors (see Ref. sothmann2015, for a review). The vast majority of the papers dealing with QD-based heat engines consider a single degenerate energy level or two non-degenerate levels weymann2013 ; donsa2014 ; sothmann2012 ; liu2011 ; kuo2013 ; wierzbicki2010 ; rossello2017 ; perroni2016 ; szukiewicz2016 ; Mazza2014 ; kuo2016 ; karol2006 ; murphy2008 ; esposito2009 ; kuo2009 ; swirkowicz2009 ; dubi2009 ; imry2010 ; leijnse2010 ; nakpathomkun2010 ; buttiker2011 ; karlstrom2011 ; trocha2012 ; muralindharan2012 ; buttiker2013 ; tseng2015 ; taylor2015 ; de2016 . The case of QDs with many levels has been addressed only in few papers agarwal2014 ; zianni2010 ; liu2010 ; kuo2010 . Moreover, the performance of QD-based heat engine has been mostly studied within the linear response regime weymann2013 ; donsa2014 ; liu2011 ; rossello2017 ; kuo2013 ; wierzbicki2010 ; Mazza2014 ; kuo2016 ; karol2006 ; murphy2008 ; kuo2009 ; dubi2009 ; imry2010 ; karlstrom2011 ; trocha2012 ; tseng2015 ; taylor2015 ; zianni2010 ; liu2010 , where the thermoelectric performance of the system is entirely characterized by $ZT$ . The case of an interacting multi-level QD beyond the linear response has not been addressed so far. On the one hand, the presence of many levels is expected to yield important consequences. Indeed, already in the linear-response regime they have an impact on the thermopower, the thermal conductance and $ZT$ bib:beenakker2 ; liu2010 ; bib:zianni_2007 ; bib:zianni_2008 . On the other, non-linear effects, relevant when larger temperature and voltage biases are applied, are important as far as power and efficiency are concerned. We emphasize that a number of experiments assessing the thermoelectric properties of QDs have been reported in Refs. llaguno2003, ; thierschmann2013, ; scheibner2007, ; godijn1999, ; reddy2007, ; dzurak1998, ; svilansPhysE2016, ; svensson2012, ; svilans2016, ; scheibner2005, ; dzurak1993, ; harman2002, ; molenkamp1992, ; staring1993, ; pogosov2006, ; scheibner2008, ; svensson2013, ; matthews2014, ; glattli2015, ; hartmann2015, ; molenkamp2015, .

In this paper we fill this gap by studying the thermoelectric properties and heat-to-work conversion performance of a multi-level QD in a multi-terminal configuration within the Coulomb blockade regime. We consider the limit of small tunnelling rates (sequential tunnelling regime) and we study both the linear and non-linear response regimes. Coulomb interaction among electrons is accounted for by a finite and small capacitance $C$ whose associated energy scale is its charging energy $(Ne)^{2}/2C$ , where $N$ is the number of electrons in the QD and $e$ is the electron charge. Moreover, we consider a generic multi-terminal structure, whereby the QD is connected to many (two or more) reservoirs. We will concentrate only on the optimization of the thermoelectric properties of the electronic system, neglecting the parasitic phononic contribution to heat transport. Our results therefore set an upper bound to the thermoelectric efficiency of the QD, approachable only in the limit in which suitable strategies to strongly reduce phonon transport are implemented.

Generalizing Refs. bib:beenakker1, and bib:beenakker2, to the multi-terminal case, by solving a set of kinetic equations one can determine the probability of occupation of the energy levels of the QD in a multi-terminal setup, thus allowing to calculate the charge and heat currents for given values of the electrochemical potentials and temperatures of the reservoirs.

In the linear response regime, where voltage and temperature biases are small, we derive closed-form expressions for the charge and heat currents and specify their limits of validity. We define local and non-local transport coefficients and express them in terms of a generating function under the assumption that the tunneling rates are energy-independent. We then derive, in the low temperature limit, analytical expressions for all transport coefficients as a function of the electrochemical potential $\mu$ . Along with the main features of the transport coefficients (located around values of $\mu$ equal to the dominant transition energies required to add or remove an electron from the QD), such expressions also describe a fine structure arising from the interplay between interaction and level quantization (controlled by the two energy scales: charging energy and level spacing). Furthermore, for the calculation of the thermal conductance we find that it is crucial to consider the presence of many levels. Within the linear response, we consider both the two-terminal and the three-terminal system aiming at addressing the performance of heat-to-work conversion. In the former case we obtain analytical expressions for the power factor $Q$ and the figure of merit $ZT$ finding, remarkably, that those quantities are simultaneously maximized for values of the electrochemical potential which differ by about $2.40k_{B}T$ with respect to the dominant transition energies. In addition, $ZT$ shows a fine structure of secondary peaks whose height is independent of the system’s parameters and can take values as large as $ZT=9$ . We compare the figure of merit with a non-interacting system, finding that Coulomb interactions dramatically increase $ZT$ by strongly suppressing the thermal conductance. For the case of three terminals with energy-independent tunneling rates, we derive analytic expressions for the maximum power and corresponding efficiency, finding that the addition of a third terminal at an intermediate temperature decreases the efficiency at maximum power, but can increase the power. We also find particular intermediate temperatures where the third terminal increases the maximum power and achieves the same efficiency of a two terminal system.

We analyze the regime beyond the linear response by numerically solving the kinetic equations, focusing on the two-terminal setup. Going beyond the linear response, i.e. considering large temperature and voltage biases, $\Delta T$ and $V$ , is interesting for various reasons. On one hand, it allows to increase the Carnot efficiency $\eta_{C}$ and the power generated by a heat engine (in our case the peak power scales approximately as $(\Delta T)^{2}$ also beyond the linear response regime). On the other hand, the efficiency at maximum power is not bounded by $\eta_{C}/2$ , as in the linear response, and can even go above the Curzon-Alhborn efficiency bib:yvon ; bib:chambadal ; bib:novikov ; bib:curzon ; bib:broeck . In literature, the scattering theory of nonlinear thermoelectric transport in quantum conductors has been developed only recently meair2013 ; whitney2013 ; Sanchez2013 . The regime beyond linear response for QD-based heat engines has been theoretically addressed in Refs. perroni2016, ; de2016, ; agarwal2014, ; buttiker2013, ; sothmann2012, ; buttiker2011, ; kuo2010, ; nakpathomkun2010, ; swirkowicz2009, ; esposito2009, ; leijnse2010, ; muralindharan2012, ; szukiewicz2016, , but limited to single or double level quantum dots.

In discussing the results, we first focus on the behavior of the non-linear Seebeck and Peltier coefficients aiming at assessing the interplay between charging energy and level spacing on these two quantities and how the Onsager reciprocity relation that connects them is violated beyond linear response. Second, we study the efficiency and output power of a heat engine. In particular, we calculate the maximum efficiency and maximum power by maximizing such quantities with respect to the applied bias voltage, for fixed values of temperature bias and electrochemical potential. The maximum efficiency shows only quantitative changes, with respect to the linear response, by increasing the temperature bias. The efficiency at maximum power instead develops peaks which go beyond the $\eta_{C}/2$ linear-response limit and approach $\eta_{C}$ for large temperature differences. Remarkably, the maximum power, normalized to its peak value, only slightly depends on the temperature bias and can be well approximated by the analytic expression obtained for the linear response regime. Moreover, we find that efficiency at maximum power and maximum power take approximately their peak values simultaneously, under the same conditions found for the linear response. Finally, we assess the impact of interactions by comparing the efficiency at maximum power in two situations, namely the case of a doubly degenerate level with interaction and the case of two non-degenerate levels without interaction. We find that, especially when charging energy and level spacing are of the order of the thermal energy, the efficiency at maximum power is much higher in the interacting case and goes above the Curzon-Alhborn efficiency. We also find that in the doubly degenerate interacting case the maximum power is enhanced by almost a factor 2 with respect to the non-degenerate case.

The paper is organized as follows: in Sec. II we describe the system under investigation and we detail the theoretical model. In Sec. III the model is specified to the linear response regime. Here we report analytic expressions obtained in the quantum limit, Sec. III.1, and calculate power and efficiency in the cases of a two-terminal system, Sec. III.2 (with the efficiencies of interacting and non-interacting QDs compared in Sec. III.3), and of a three-terminal system, Sec. III.4. Sec. IV is devoted to the discussion of the regime beyond the linear response in a two-terminal system: in Sec. IV.1 we study the non-linear Seebeck and Peltier coefficients and in Sec. IV.2 we discuss efficiency and output power of a heat engine. Finally, in Sec. V we draw our conclusions and discuss future developments. In addition, the paper includes four appendices where the details of some calculations and the analytic non linear study of a single energy level QD system are reported.

II Multi-level interacting QD

In this Section we briefly outline the formalism used to describe the thermoelectric properties of a multilevel interacting QD. We will only consider electron transport, neglecting any contribution due to phonons. As shown in Fig. 1 (top), the QD is tunnel-coupled to ${\cal N}$ electron reservoirs, each characterized by a given temperature $T_{\alpha}$ and electrochemical potential $\mu_{\alpha}$ , so that the occupation of the electrons within reservoir $\alpha$ follows the Fermi distribution

[TABLE]

where $k_{B}$ is Boltzmann’s constant. In Fig. 1 (bottom), $E_{p}$ (with $p=1,2,\dots$ labeled in ascending order) are the QD single-electron energy levels. These levels can be shifted by means of an applied gate voltage.

The QD is weakly coupled to the reservoirs through large tunneling barriers. More precisely, we assume that thermal energy $k_{B}T$ , level spacing and charging energy are much larger than the coupling energy between reservoirs and QD [ $\hbar\sum_{\alpha}\Gamma_{\alpha}(p)$ , where $\Gamma_{\alpha}(p)$ is the tunneling rate from level $p$ to reservoir $\alpha$ , which we assume independent of the number $N$ of electrons inside the dot]. As a consequence, the charge on the QD is quantized, i.e. each energy level $E_{p}$ can have either zero or one electron, $n_{p}=0$ or $n_{p}=1$ (any degeneracy, like electron spin, can be taken into account counting each level multiple times), and transport occurs due to single-electron tunneling processes (sequential tunneling regime). The electrostatic energy associated with the electrons within the QD is given by $U(N)=E_{C}N^{2}$ , where $E_{C}=e^{2}/2C$ , $N=\sum_{i}n_{i}$ is the total number of electrons within the QD, and $C$ is the capacitance of the QD. The QD is described by states characterized by a set of occupation numbers $\{n_{i}\}$ relative to the energy levels. The QD changes state whenever a single-electron tunneling process takes place. The non-equilibrium probability for a given state $\{n_{i}\}$ to occur, $P(\{n_{i}\})$ , can be computed bib:beenakker1 ; bib:nazarov by writing a straightforward set of balance equations for $P\left(\{n_{i}\}\right)$ . Our aim is to compute in stationary conditions the charge, energy and heat currents out of the electron reservoirs [denoted in Fig. 1 (top) by $J_{\alpha}^{c}$ , $J_{\alpha}^{u}$ and $J_{\alpha}^{h}$ , respectively], induced by the temperature and electrochemical potential differences.

II.1 Kinetic equations

In what follows we describe a generalization of the method put forward by Beenakker in Refs. bib:beenakker1, and bib:beenakker2, . The single-electron tunneling processes that contribute to changing over time the probability $P(\{n_{i}\})$ are due to electrons that tunnel from the QD to the reservoirs and vice-versa. For an electron exiting the QD, initially with $N$ electrons, from energy level $E_{p}$ and going into reservoir $\alpha$ at energy $E^{\text{fin}}$ , energy conservation imposes that

[TABLE]

On the contrary, for an electron that tunnels from an initial state in reservoir $\alpha$ at energy $E^{\text{in}}$ to the level $E_{p}$ in the QD that initially had $N$ electrons, energy conservation imposes that

[TABLE]

$P(\{n_{i}\})$ can then be determined by the following set of kinetic equations, one for each configuration $\{n_{i}\}$ :

[TABLE]

where we have introduced the notation

[TABLE]

and

[TABLE]

for the QD states. The first term in Eq. (4) accounts for the decrease of the probability $P(\{n_{i}\})$ , with the QD initially in the state $\{n_{i}\}$ , due to an electron coming from a reservoir and occupying an empty level in the QD. The rate of electrons coming from reservoir $\alpha$ will be given by a sum over all empty levels $p$ (such that $n_{p}=0$ ) of the tunnel rate $\Gamma_{\alpha}(p)$ , multiplied by the probability of finding the QD in this state, $P\left(\{n_{i}\}\right)$ , and multiplied by the reservoir’s occupation $f_{\alpha}\left(E^{\text{in}}(N)\right)$ at the correct energy $E^{\text{in}}(N)$ to tunnel to level $p$ . The second term accounts for the decrease of the probability $P(\{n_{i}\})$ , with the QD initially in the state $\{n_{i}\}$ , due an electron leaving the QD from an occupied level to tunnel into a reservoir. The third term accounts for the increase of the probability $P(\{n_{i}\})$ if the QD is in a state with an extra electron in level $p$ with respect to $\{n_{i}\}$ , and if this electron leaves the QD, tunneling to the reservoirs. The forth term accounts for the increase of the probability $P(\{n_{i}\})$ if the QD is in a state with a missing electron in level $p$ with respect to $\{n_{i}\}$ , and if this electron enters the QD in level $p$ , tunneling from the reservoirs. The stationary solution of the kinetic equations, obtained imposing $\partial P/\partial t=0$ , together with the normalization request

[TABLE]

provides a complete set of equations that uniquely defines $P\left(\{n_{i}\}\right)$ . The sum over $\{n_{i}\}$ means the sum over $n_{i}=0,1$ , with $i=1,2,...$ .

II.2 Charge, energy, and heat currents

Charge $J^{c}_{\alpha}$ and energy $J^{u}_{\alpha}$ currents flowing from reservoir $\alpha$ to the QD can be calculated as the sum of all possible tunneling processes, since the QD can be in any state $\{n_{i}\}$ with probability $P(\{n_{i}\})$ and an electron can tunnel into or out of any energy level $E_{p}$ . More precisely, for the charge current we have

[TABLE]

$e$ being the electronic charge, while for the energy current we have

[TABLE]

$E^{\text{in}}(N)$ [ $E^{\text{fin}}(N)$ ] being the energy carried by an electron entering (exiting) the QD. The heat currents exiting the reservoirs can be calculated as $J^{h}_{\alpha}=J^{u}_{\alpha}-\frac{\mu_{\alpha}}{e}J^{c}_{\alpha}$ . Using Eqs. (8) and (9), we find that

[TABLE]

In order to numerically determine the stationary probability distribution $P(\{n_{i}\})$ from the kinetic equations, we will consider a finite number $L$ of energy levels. 111The results can be interpreted as the exact solution of a system with $L$ energy levels or as an approximate solution of a system with infinite levels. In the latter case, one needs to check that the result are stable with increasing number of energy levels.. By organizing the values of $P(\{n_{i}\})$ into a $2^{L}$ -components vector $\vec{P}$ (two choices $n_{i}=0,1$ for each level), the kinetic equations (4) for the stationary probability distribution, $\frac{\partial\vec{P}}{\partial t}=\vec{0}$ , can be represented as the homogeneous linear system $M\vec{P}=\vec{0}$ , where $M$ is a $2^{L}\times 2^{L}$ matrix. $M$ must have a null space of at least dimension 1, otherwise the only possible solution would be the trivial one ( $\vec{P}\equiv 0$ ). This is demonstrated in App. A by showing that summing together all the stationary kinetic equations yields zero. We can thus find the probabilities by including the normalization condition, Eq. (7).

By defining $\widetilde{N}=\sum\limits_{i\neq p}n_{i}$ it is possible to show that the kinetic equations for the stationary probability distribution can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

To derive Eq. (11) we have used the fact that $\widetilde{N}=N$ , for the terms in the kinetic equations proportional to $\delta_{n_{p},0}$ , and $\widetilde{N}=N-1$ , for the terms in the kinetic equations proportional to $\delta_{n_{p},1}$ , and the identity

[TABLE]

stemming from Eqs. (2) and (3). It is worth mentioning that using the kinetic equations in the form of Eq. (11), it is possible to prove that the kinetic equations always allows a non-trivial solution (see App. A).

II.3 Detailed balance equations

It is clear that the kinetic equations (11) are automatically satisfied when the following set of equations

[TABLE]

is fulfilled for all values of $p$ and for all sets of occupation numbers $\{n_{i}\}$ . Following Ref. bib:beenakker1, , Eqs. (15) are hereafter referred to as detailed balance equations (DBEs). Eqs. (15) represent a set of $L\cdot 2^{L-1}$ equations, since $p$ can take $L$ values and, at a given $p$ , all other occupation numbers ( $n_{1},\dots,n_{p-1},n_{p+1},\dots,n_{L}$ ) can be chosen in $2^{L-1}$ different ways. Of course, if a solution to the DBEs exists, than it is also a solution of the kinetic equations. We can show, however, that the DBEs are not in general consistent, i.e. no set of $P(\{n_{i}\})$ exists that can simultaneously satisfy all the DBEs (see App. B). This is also true in the linear response regime. In this case, however, we could prove (see App. C) that the DBEs are consistent if $E_{C}=0$ , or if $\Delta T_{\alpha}=0$ for all $\alpha$ , or when the tunneling rates are proportional to each other, namely when $\Gamma_{\alpha}(p)=k_{\alpha}\Gamma_{1}(p)$ , for $\alpha>1$ , $k_{\alpha}$ being constants. Note that this condition is trivially satisfied when the rates $\Gamma_{\alpha}$ do not depend on $p$ . As a result, the DBEs do not allow in general a solution, but when they do they are useful in computing analytically the energy and heat currents in the linear response regime (see Sec. III).

II.4 Level balance equations

We will now derive a set of equations that is always consistent and that can be used in the general case to obtain a closed-form expression of the charge current in the linear response regime (see Sec. III). We impose that, in stationary conditions, the rate of electrons entering any given QD energy level $p$ must equal the rate of electrons leaving that energy level. For electrons tunneling into the QD, initially with $N$ electrons, from any reservoir, one has to require that level $p$ is empty and must consider all possible states $(\{n_{i}\},n_{p}=0)$ , where $\widetilde{N}=N$ . The total rate of electrons entering energy level $E_{p}$ is given by

[TABLE]

For electrons tunneling out of the QD, initially with $N$ electrons, to any reservoir one has to require that level $p$ is occupied and must consider states with $(\{n_{i}\},n_{p}=1)$ , where $\widetilde{N}=N-1$ . The total rate of electrons leaving level $E_{p}$ is given by

[TABLE]

If we equate the rates of electron entering and leaving level $p$ , and use the notation introduced in Eqs. (12) and (13), we obtain

[TABLE]

We will refer to this set of $L$ equations (one for each energy level $p$ ) as the level balance equations (LBEs) 222Note that the DBEs (15) could be obtained using the same arguments, but assuming that the QD is in a specific configuration $(\{n_{i}\})_{i\neq p}$ , instead of summing over all possible such configurations. To obtain the LBEs, the probability that level $p$ is occupied is given by the marginal probability $\sum\nolimits_{\{n_{i}\}_{i\neq p}}P\left(\{n_{i}\},n_{p}=1\right)$ .. Note that, using an argument similar to that put forward in App. A, it is possible to prove that Eqs. (18) can be obtained from the kinetic equations, thus the LBEs are always consistent with the kinetic equations. However, the number of LBEs (equal to $L$ ) is smaller than the number of kinetic equations (equal to $2^{L}$ ) and they might not be sufficient to determine the probabilities $P\left(\{n_{i}\}\right)$ . Note that one can prove that Eqs. (18) yield charge current conservation: $\sum_{\alpha}J^{c}_{\alpha}=0$ .

II.5 Output power and efficiency

Under steady-state conditions, the output power $P$ of a multi-terminal system is given by the sum of all the heat currents

[TABLE]

If $P>0$ , the system behaves as a heat engine, i.e. converting heat into work. In this situation the efficiency $\eta$ is defined as the ratio between the output power and the total heat current absorbed by the system

[TABLE]

where the sum over $\alpha^{\prime}$ runs over all positive heat currents. For a two-terminal system the efficiency cannot exceed the Carnot efficiency defined as $\eta_{C}=1-T_{1}/T_{2}$ , with $T_{2}>T_{1}$ . In addition, for a multi-terminal system $\eta$ cannot go beyond the two-terminal Carnot efficiency Mazza2014 calculated using the hottest and coldest temperatures among $T_{1},T_{2},...,T_{\cal N}$ .

We define the temperature and electrochemical potential differences as $T_{\alpha}=T+\Delta T_{\alpha}$ and $\mu_{\alpha}=\mu+\Delta\mu_{\alpha}$ , with $\alpha=1,...,{\cal N}$ , and choosing reservoir $1$ as the reference value, i.e. $\Delta T_{1}=\Delta\mu_{1}=0$ . In what follows we fix the values of $\Delta T_{\alpha}$ and calculate the maximum output power, $P_{\text{max}}$ , and maximum efficiency, $\eta_{\text{max}}$ , by varying $\Delta\mu_{\alpha}$ . We also consider the efficiency at maximum power, $\eta(P_{\text{max}})$ , which is the efficiency when the values of $\Delta\mu_{\alpha}$ are chosen to maximize the power.

For a two-terminal system within the linear response regime, i.e. when the charge and heat currents depend linearly on the temperature and electrochemical potential differences, both the output power and efficiency can be written in terms of the transport coefficients, namely the electrical conductance $G$ , the thermopower $S$ and the thermal conductance $K$ , which will be defined in Eqs. (34), (35) and (36) by setting $\alpha=\beta=2$ . Defining $\Delta T\equiv\Delta T_{2}>0$ , we have the following relations PRreview ; bib:ioffe2

[TABLE]

where $Q=GS^{2}$ is the power factor and $ZT=GS^{2}T/K$ is the (dimensionless) figure of merit. As we can see in Eqs. (22) and (23), both $\eta(P_{\text{max}})$ and $\eta_{\rm max}$ are monotonous growing functions of $ZT$ ; the only restriction imposed by thermodynamics is $ZT\geq 0$ . When $ZT=0$ both $\eta_{\rm max}$ and $\eta(P_{\text{max}})$ vanish, while for $ZT\to\infty$ , $\eta_{\rm max}\to\eta_{C}$ , and $\eta(P_{\text{max}})\to\eta_{CA}$ , where $\eta_{CA}=\eta_{C}/2$ is the so-called Curzon-Ahlborn efficiency bib:yvon ; bib:chambadal ; bib:novikov ; bib:curzon ; bib:broeck in linear response.

III Linear response regime

As already mentioned above, in the linear response regime the applied temperature and electrochemical potential biases are small enough so that the currents depend linearly on them. Assuming that $|\Delta T_{\alpha}|\ll T$ and $|\Delta\mu_{\alpha}|\ll k_{B}T$ , we follow Refs. bib:beenakker1, and bib:beenakker2, and suppose the probability $P(\{n_{i}\})$ to differ from its equilibrium distribution $P_{\text{eq}}(\{n_{i}\})$ in the following way:

[TABLE]

where $\psi$ is a “small” function. In Eq. (24)

[TABLE]

is the Gibbs distribution in the grand canonical ensemble, when all reservoirs have the same temperature and electrochemical potential, with grand partition function given by

[TABLE]

In our expressions we will consider terms up to first order in $\psi$ , $\Delta T_{\alpha}/T$ , and $\Delta\mu_{\alpha}/k_{B}T$ . By linearizing the LBEs with respect to the above small quantities, Eq. (18), one finds the relation

[TABLE]

where $f(E)$ stands for the Fermi distribution at temperature $T$ and electrochemical potential $\mu$ . By expressing $P(\{n_{i}\})$ in terms of $\psi(\{n_{i}\})$ and linearizing Eq. (8), we can use Eq. (27) to remove $\psi(\{n_{i}\})$ from the charge current, and we find the following closed-form expression:

[TABLE]

where $\Gamma_{\text{tot}}(p)=\sum_{\alpha}\Gamma_{\alpha}(p)$ , and

[TABLE]

is the energy needed to add to the QD, initially occupied by $N-1$ electrons, the $N$ -th electron to level $p$ (and equivalently for the inverse process). In Eq. (28) the quantity

[TABLE]

is the equilibrium probability of finding $N$ electrons in the QD, and

[TABLE]

is the equilibrium conditional probability of finding level $p$ occupied, when $N$ electrons are in the QD. Note that expression (28) coincides with the one that can be derived using the DBEs. However, the above derivation which uses the LBEs shows that Eq. (28) is always valid within the linear response regime.

Unfortunately, we were able to derive a closed-form expression for the energy current using Eq. (27) only in the absence of interaction ( $E_{C}=0$ ). For $E_{C}\neq 0$ , the energy current $J_{\alpha}^{u}$ can be derived using the relation

[TABLE]

obtained by linearizing the DBEs, Eq. (15) (which is equivalent to removing the sum over $\{n_{i}\}$ from Eq. (27)). Thus, in the domain of validity of the DBEs, the heat current can be written as

[TABLE]

We can now define the transport coefficients, namely the electrical conductance $G_{\alpha\beta}$ , the thermopower $S_{\alpha\beta}$ and the thermal conductance $K_{\alpha\beta}$ for the multi-terminal case, as Mazza2014

[TABLE]

and

[TABLE]

Note that index $\beta$ takes values in the range $2,...,{\cal N}$ , since reservoir 1 is chosen as the reference. Local and non-local transport coefficients are distinguished depending on whether the two indices are, respectively, equal or different.

The expressions for the currents [(28) and (33)] have an intuitive interpretation. Indeed, the currents depend on the probability that a given energy level of the QD is occupied [ $P_{\text{eq}}(N)F_{\text{eq}}(E_{p}|N)$ ] times the probability that there is an empty state with the correct energy in the reservoir [ $1-f(\epsilon(N,p))$ ]. The sum over all energy levels $p$ and over the total number of electrons $N$ in the QD accounts for all the various tunneling processes that can occur. Moreover, as far as the energy current is concerned, $\epsilon(N,p)$ is the energy carried by an electron that leaves the QD from level $p$ when $N$ electrons are present before the tunneling process, or equivalently, $\epsilon(N,p)$ is the energy carried by an electron that enters the QD into level $p$ increasing the number of total electrons to $N$ . We recall that Eq. (28) for the charge current is always valid, while Eq. (33) holds only when the DBEs are valid (see App. C). The expressions for the charge and heat currents, in the case of two terminals, coincide with the ones obtained in Refs. bib:beenakker1, ; bib:beenakker2, ; bib:zianni_2007, .

If we assume that the tunneling rates do not depend on $p$ , i.e. $\Gamma_{\alpha}(p)=\Gamma_{\alpha}$ , we can rewrite the charge and heat currents in Eqs. (28) and (33) as follows:

[TABLE]

and

[TABLE]

where ${\cal P}$ is the linear functional

[TABLE]

with

[TABLE]

For a two terminal system it is possible to define analogous equations that do not require the tunneling rates to be energy-independent paolo2016 , but we will not consider this case. We can thus use the definitions of the transport coefficients given in Eqs. (34), (35) and (36) to write them in terms of the functional ${\cal P}$ as

[TABLE]

These expressions only require the calculation of ${\cal P}[(\epsilon-\mu)^{k}]$ , with $k=0,1,2$ , and make manifest various properties of the transport coefficients. Namely, i) all three transport coefficients are symmetric matrices (as required by the Onsager relations in the presence of time-reversal symmetry); ii) $G_{\alpha\beta}$ and $K_{\alpha\beta}$ have non-local terms, while the thermopower is only local (non-zero non-local $S_{\alpha\beta}$ occur when relaxing the assumption for which the tunneling rates do not depend on the energy levels Mazza2014 ); iii) $\sum\nolimits_{\alpha}G_{\alpha\beta}=\sum\nolimits_{\alpha}K_{\alpha\beta}=0$ , stemming from charge and energy conservation.

By defining the generating function

[TABLE]

we can write the transport coefficients as follows:

[TABLE]

In the next Subsection, we will compute an analytic expression for $\Omega[\lambda]$ in the quantum limit.

III.1 Quantum limit

The quantum limit is characterized by having the energy spacing between levels of the QD and the charging energy much bigger than $k_{B}T$ [while $k_{B}T\gg\hbar\Gamma_{\alpha}(p)$ ]. We start by observing that the sum over $p$ and $N$ in Eq. (39) accounts for the fact that electrons can enter or leave the QD with energy $\epsilon(N,p)$ through, in principle, any energy level $E_{p}$ with the QD being occupied by any number of electrons $N$ . The transition energies $\epsilon(N,p)$ are schematically shown in Fig. 2.

At low temperatures we expect the lowest energy levels of the QD to be occupied, so that, if there are initially $N-1$ electrons in the QD, electrons can flow mainly through level $p=N$ . Such process gives the dominant contribution to transport and is represented by the dominant transition energy $\tilde{\epsilon}(N)\equiv\epsilon(N,p=N)$ (depicted in bold in Fig. 2). Therefore, in the quantum limit one expects to get the sum over $N$ and $p$ appearing in ${\cal P}$ [Eq. (39)] reduced to few dominant terms [the three equilibrium probabilities $P_{\text{eq}}(N)$ , $F_{\text{eq}}(E_{p}|N)$ and $f(\epsilon(N,p))$ becoming very sharp functions].

Following Ref. bib:beenakker2, , one finds (see Ref. paolo2016, for details) that the dominant contribution to $P_{\text{tot}}$ in Eq. (39) occurs when $N=N_{\text{min}}$ is the integer that minimizes the quantity

[TABLE]

and for values of $p$ such that $\epsilon(N_{\text{min}},p)$ is between the electrochemical potential $\mu$ and the dominant transition energy, i.e. such that $\tilde{\epsilon}(N_{\text{min}})\leq\epsilon(N_{\text{min}},p)<\mu$ or $\mu<\epsilon(N_{\text{min}},p)\leq\tilde{\epsilon}(N_{\text{min}})$ . In the former case, $p=N_{\text{min}},N_{\text{min}}+1,N_{\text{min}}+2,...,\bar{p}$ , where $\bar{p}$ is the largest integer such that $\epsilon(N_{\text{min}},p)<\mu$ . In the latter case, $p=N_{\text{min}},N_{\text{min}}-1,N_{\text{min}}-2,...,\bar{p}$ where $\bar{p}$ is the smallest integer such that $\epsilon(N_{\text{min}},p)>\mu$ . We then find that Eq. (39) becomes

[TABLE]

where we have defined $\Delta_{\text{min}}\equiv\tilde{\epsilon}(N_{\text{min}})-\mu$ . Eq. (45) only keeps the dominant terms in the low-temperature limit.

This approximation must be improved when $\bar{p}=N_{\text{min}}$ , that is, the sum in ${\cal P}$ in Eq. (45) reduces to the single term $p=N_{\text{min}}$ and transport is provided by the dominant transition energy $\tilde{\epsilon}(N_{\text{min}})$ only. If in this case one imposes that $J^{c}_{\alpha}=0$ , then one obtains $J^{h}_{\alpha}=0$ , since $J^{h}_{\alpha}\propto J^{c}_{\alpha}$ . As a consequence one gets $K_{\alpha\beta}=0$ 333This approximation for two terminals leads to zero thermal conductance and to the Carnot efficiency, as discussed for single-level quantum dots bib:mahan ; bib:linke ; bib:humphrey .. Thus, when $\bar{p}=N_{\text{min}}$ , we improve our approximation of ${\cal P}$ by extending the sum over $p$ to the two nearest integers, $p=N_{\text{min}}\pm 1$ . We have numerically verified that this approximation is valid when $2E_{C}>\Delta E$ .

In order to obtain analytical expressions for the transport coefficients, hereafter we focus on the case of equidistant levels, $E_{p}-E_{p-1}=\Delta E$ . After introducing the parameter

[TABLE]

and defining the integer $N_{J}\equiv\bar{p}-N_{\text{min}}$ , we obtain

[TABLE]

Using Eqs. (43) and (47), we finally obtain the multiterminal transport coefficients

[TABLE]

We have computed $G_{\alpha\beta}$ and $S_{\alpha\beta}$ setting $\xi=0$ , since the term proportional to $\xi$ in Eq. (47) only yields minor corrections that make these quantities more “smooth” as a function of $\mu$ ; the calculation of $K_{\alpha\beta}$ instead requires a non null value of $\xi$ . Eqs. (48) exhibit a number of interesting features. First of all, $G_{\alpha\beta}$ shows peaks as a function of the electrochemical potential $\mu$ every time $\Delta_{\text{min}}=0$ , namely when $\mu=\tilde{\epsilon}(N)$ . For example, the $N$ -th peak corresponds to $\mu$ equal to the $N$ -th dominant transition energy,

[TABLE]

We set $E_{p}=(p-1)\Delta E$ , for $p=1,2,\dots$ , and therefore the separation between two nearby peaks is given by $\Delta E+2E_{C}$ . Due to the factor $\cosh^{-2}(\frac{\Delta_{\text{min}}}{2k_{B}T})$ in $G_{\alpha\beta}$ , these peaks have a bell shape with amplitude of the order of $k_{B}T$ . On the other hand, the thermal conductance $K_{\alpha\beta}$ has plateaus of width $2\Delta E$ around $\mu_{N}$ bib:zianni_2007 , corresponding to the second line of the expression of $K_{\alpha\beta}$ for $N_{J}=0$ . For $N_{J}\neq 0$ , the thermal conductance $K_{\alpha\beta}$ is then exponentially suppressed due to the term $e^{-|\Delta_{\text{min}}|/k_{B}T}$ . The local thermopower $S_{\alpha\alpha}$ vanishes at the values $\mu_{N}$ where the electrical and thermal conductances $G_{\alpha\alpha}$ and $K_{\alpha\alpha}$ exhibit a maximum. $S_{\alpha\alpha}$ has a linear dependence on $\mu$ with slope $dS_{\alpha\alpha}/d\mu=-1/eT$ , with jumps when either $N_{\rm min}$ or $N_{J}$ change by one. Therefore we have, for $E_{C}\gg\Delta E$ , main oscillations of period $\Delta E+2E_{C}$ , and a fine structure with spacing $\Delta E$ bib:beenakker2 . We note that the fine structure is present also for $G_{\alpha\beta}$ and $K_{\alpha\beta}$ , but in these cases the amplitude of the fine structure oscillations is exponentially small.

III.2 Two-terminal system

In the two-terminal case, the matrices $G_{\alpha\beta}$ , $S_{\alpha\beta}$ , and $K_{\alpha\beta}$ reduce to the familiar transport coefficients, namely the electrical conductance $G=G_{22}$ , the thermopower $S=S_{22}$ , and the thermal conductance $K=K_{22}$ . From Eqs. (48) we recover the formulas for $G$ and $S$ well-known in literature bib:beenakker2 , while our expression for $K$ coincides for $N_{J}=0$ with the result of Ref. bib:zianni_2007, , but also provides the fine structure oscillations for $N_{J}\neq 0$ . Although such oscillations are not appreciable in $K$ as a function of $\mu$ , we will see below that they give rise to a visible fine structure in $ZT$ .

As we have shown in Sec. II.5, within the linear response regime the relevant quantities to characterize the performance of a thermoelectric device are the power factor $Q$ and the figure of merit $ZT$ . From the expressions of the transport coefficients in Eqs. (48), specified for the two-terminal case, one can compute $Q$ and $ZT$ analytically within the quantum limit. The obtained expressions are given below, and compared with a numerical calculation performed using the kinetic equations.

III.2.1 Power factor

Let us start by studying the power factor $Q$ . Within the quantum limit, we find that

[TABLE]

where we have defined $\gamma\equiv\Gamma_{1}\Gamma_{2}/(\Gamma_{1}+\Gamma_{2})$ .

As for $G$ , the power factor $Q$ is dominated by a fast decrease, given by the term $\cosh^{-2}{(\Delta_{\text{min}}/2k_{B}T)}$ , thus becoming vanishingly small within a few $k_{B}T$ around $\mu=\mu_{N}$ [see Fig. 3(a)]. In fact the fine structure, given by the terms with $N_{J}$ , is not visible in Fig. 3(a) due to the rapid suppression given by the $\cosh^{-2}{(\Delta_{\text{min}}/2k_{B}T)}$ term in Eq. (50). Differently from $G$ , the power factor vanishes at $\mu=\mu_{N}$ , due to the fact that in this point the thermopower $S=0$ . So as $\mu$ moves away from $\mu_{N}$ , $Q$ increases quadratically due to the linear growth of the thermopower with $\mu$ , and then it rapidly decreases within a few $k_{B}T$ due to the $\cosh^{-2}(\Delta_{\text{min}}/2k_{B}T)$ term. Hence there are two symmetric peaks around $\mu=\mu_{N}$ , within a few $k_{B}T$ . These double peaks are the dominant feature of Fig. 3(a) and identify the optimal values of $\Delta_{\text{min}}$ (and consequently of $\mu$ ) to obtain the absolute maximum power $P_{\text{peak}}$ , namely when the power factor $Q$ is maximum, $Q=Q^{*}$ . From Eq. (50), we obtain that $Q$ is maximum for values $\Delta_{\text{min}}^{*}$ of $\Delta_{\text{min}}$ such that

[TABLE]

The numerical solution is $\Delta_{\text{min}}^{*}\simeq\pm 2.40k_{B}T$ , which corresponds to $\mu=\mu_{N}\pm 2.40k_{B}T$ . This result does not depend on any energy scale of the system except for $k_{B}T$ , and coincides with the non-interacting single-level case (see App. D). The value $Q^{*}$ of $Q$ in these points is

[TABLE]

so that the peaks of the power factor only depend on $\gamma$ and on the reference temperature. In conclusion, if we want to extract maximum power from this system, we must choose $\mu=\mu_{N}\pm 2.40k_{B}T$ . We will now show that also $ZT$ reaches a maximum at these same values of the electrochemical potential confirming that these are the optimal values for heat to work conversion in the quantum limit linear response regime.

III.2.2 Figure of merit

Let us now study the figure of merit $ZT$ in the quantum limit. To obtain a more manageable analytical expression, we compute $K$ from the function $\Omega[\lambda]$ expanded to the first order in $\xi$ [this corresponds to approximating $K$ with a constant plateau when $N_{J}=0$ , namely $4\cosh^{2}(\Delta_{\text{min}}/2k_{B}T)$ is neglected with respect to $e^{\Delta E/k_{B}T}$ in the last line of Eqs. (48)]. We then obtain

[TABLE]

which implies that the behavior of $ZT$ is different for the two cases $|\mu-\mu_{N}|<\Delta E$ ( $N_{J}=0$ ) and $|\mu-\mu_{N}|>\Delta E$ ( $N_{J}\neq 0$ ). In the former case $K$ exhibits a plateau, so that $ZT$ is directly proportional to $Q$ and therefore it has the same double peak structure at $\mu=\mu_{N}\pm 2.40k_{B}T$ . This is clearly shown in Fig. 3(b) where $ZT$ is plotted as a function of $\mu$ .

The value of $ZT$ in these points is

[TABLE]

This results has been obtained also in Ref. Tsaousidou2010, . Eq. (54) shows that in the limit $\Delta E/k_{B}T\rightarrow\infty$ , we have that $ZT\rightarrow\infty$ . For example, for $\Delta E=6k_{B}T$ , we reach $ZT^{*}\approx 5$ ; for $\Delta E=10k_{B}T$ , we reach $ZT^{*}\approx 97$ , and so on. This is consistent with Mahan and Sofo’s observation bib:mahan that a narrow transmission function yields $ZT\rightarrow\infty$ . Furthermore, these peaks in $ZT$ correspond to peaks in $Q$ , so in these points we can maximize $P_{\text{max}}$ and $\eta(P_{\text{max}})$ simultaneously. Instead, when $N_{J}\neq 0$ , $ZT$ has a discontinuity every time $\mu=\epsilon(N,p)$ with $p\neq N$ , which means with a $\Delta E$ spacing. This fine structure is the origin of the saw-tooth oscillations of Fig. 3(b). The value of $ZT$ in each $\mu=\epsilon(N,p)$ is given by

[TABLE]

The height of these peaks, as opposed to $ZT^{*}$ , has no dependance on the parameters of the system. The highest peak is obtained for $|N-p|=2$ , where $ZT_{p=N\pm 2}=9.$ For values of $p$ distant from $N$ , the height of the peak decreases to an asymptotic value of $ZT_{\infty}=3$ .

III.3 Comparison with a non-interacting QD

Here we compare the efficiency of an interacting QD (with $2E_{C}>\Delta E$ ) with the efficiency of a non-interacting QD ( $E_{C}=0$ ) that has the same energy spacing $\Delta E$ ; the comparison is performed within the linear response quantum limit for a two-terminal setup. The generating function in Eq. (47) cannot be used in the case $E_{C}=0$ since it requires $2E_{C}>\Delta E\gg k_{B}T$ . The generating function for $E_{C}=0$ will be denoted as $\Omega_{\text{NI}}[\lambda]$ (where NI stands for “non interacting”) and calculated as follows. As we can see from Eq. (42), we must compute ${\cal P}[\exp{\{\lambda(\epsilon(N,p)-\mu)\}}]$ using the definition of ${\cal P}$ given in Eq. (39). When $E_{C}=0$ , $\epsilon(N,p)=E_{p}$ , so the transition energies correspond to the energy levels of the QD, and they do not depend on the number of electrons in the QD. As a consequence, there is no dependance on $N$ in the argument of ${\cal P}$ , so we can explicitly perform the sum over $N$ in Eq. (39) yielding the following expression:

[TABLE]

Each term in the series in Eq. (56), as a function of $\mu$ , is a bell-shaped function centered around $E_{p}$ of width set by $k_{B}T$ . Therefore, within the quantum limit we can restrict the sum over $p$ to the three bell-shaped functions closest to $\mu$ , namely $p=N_{\text{min}},N_{\text{min}}\pm 1$ ( $N_{\text{min}}$ , defined as in Eq. (44), is such that $E_{N_{\text{min}}}$ is the closest energy level to $\mu$ ). This approximation allows us to find

[TABLE]

where $\Delta_{\text{min}}$ in the non-interacting case reduces to the distance between $\mu$ and the nearest energy level: $\Delta_{\text{min}}=E_{N_{\text{min}}}-\mu$ .

Let us now compare Eq. (57) with its interacting counterpart, Eq. (47), setting $N_{J}=0$ , thus neglecting the fine structure. We notice that the two expressions are identical when $\xi=0$ ; this implies that $G_{\alpha\beta}$ and $S_{\alpha\beta}$ are equal in the interacting and non-interacting case [in Eqs. (48) $G_{\alpha\beta}$ and $S_{\alpha\beta}$ are calculated at $\xi=0$ ], while $K_{\alpha\beta}$ is different in the two cases, since it is determined by the term proportional to $\xi$ in Eqs. (47) and (57). We find that

[TABLE]

A comparison between the interacting and non-interacting thermal conductances is plotted for a two-terminal system in Fig. 4, using Eqs. (48) and (58) at equal $\Delta E=10k_{B}T$ .

In this figure the dominant transition energy of the interacting system is located, see Eq. (49), at $\mu=\mu_{N=2}=160\,k_{B}T$ , while the energy levels of the non-interacting system are chosen as $E_{p}=(p-1)10k_{B}T$ so that in both cases $G$ has a peak in $\mu=160\,k_{B}T$ . We have verified that a numerical calculation is in very good agreement with Eq. (58) using the parameters of Fig. 4. As we can see from Fig. 4, $K^{(\text{NI})}=K^{(\text{NI})}_{22}$ and $K=K_{22}$ are very different. While $K$ has a plateau of width $2\Delta E$ centered in $\mu_{N=2}$ , $K^{(\text{NI})}$ has a minimum in $\mu_{N=2}$ , and reaches a maximum value for $\mu$ between two energy levels. Furthermore, by comparing Eq. (58) with the last line of Eqs. (48), we see that $K^{(\text{NI})}\approx 2\cosh{\left(\Delta_{\text{min}}/k_{B}T\right)}K$ , so the minimum of $K^{(\text{NI})}$ , occurring at $\Delta_{\text{min}}=0$ , is twice the maximum of $K$ , and $K^{(\text{NI})}$ increases exponentially respect to $K$ as $\Delta_{\text{min}}$ increases.

Intuitively, the striking difference between the two models can be explained as follows. As discussed in App. D, if we consider a single energy level QD the thermal conductance vanishes ( $K=0$ ) since $K$ is computed at zero charge current and charge and heat currents are proportional in this case. However, $K$ can be finite when at least two energy levels are available, and gets bigger by increasing the flux of electrons tunneling at different energies. Now, Coulomb interaction produces a correlation between electrons tunneling at different energies. Namely, if one electron enters the QD the electrostatic energy increases by $2E_{C}$ , preventing other electrons from entering the QD at any other energy level. Therefore, until that electron tunnels out of the QD, all other processes are blocked: this is a manifestation of Coulomb blockade. On the contrary, in the non-interacting case all tunneling events are independent. This correlation is thus responsible for suppressing simultaneous tunneling through different energy levels in the interacting case, which results in a suppression of $K$ .

So in general $K$ is much smaller than $K^{(\text{NI})}$ . As a consequence of these observations, we expect $ZT$ to be suppressed in the non-interacting case.

By setting $\xi=0$ in Eq. (58), we find

[TABLE]

Comparing Eq. (59) with Eq. (53), we see that, for $N_{J}=0$ , $ZT=2\cosh{\left(\Delta_{\text{min}}/k_{B}T\right)}ZT_{\text{NI}}$ , so $ZT_{\text{NI}}$ is exponentially suppressed as $\Delta_{\text{min}}$ increases. Given this suppression, the maximum of $ZT_{\text{NI}}$ occurs at $\Delta_{\text{min}}^{*}\approx\pm 1.36k_{B}T$ , corresponding to

[TABLE]

This value is approximately 6 times smaller than $ZT^{*}$ for the interacting case, see Eq. (54).

Furthermore, since $K$ is “flat” around the dominant transition energies in the interacting case, the peak power $P_{\text{peak}}$ and the figure of merit $ZT^{*}$ are reached at the same electrochemical potential, $\Delta_{\text{min}}^{*}\approx\pm 2.40k_{B}T$ , so $ZT(P_{\text{peak}})=ZT^{*}$ . Instead in the non-interacting case these two quantities are not simultaneously maximized, due to the strong dependance of $K^{(\text{NI})}$ on $\mu$ , so we have that

[TABLE]

which is approximately 11 times smaller than $ZT(P_{\text{peak}})$ , see Eq. (54).

In conclusion, within the linear response quantum limit, an interacting QD (with $2E_{C}>\Delta E$ ) has a considerably higher $ZT$ with respect to a non-interacting QD both at peak efficiency, and at peak power, while having the same power factor ( $G$ and $S$ being equal).

At last we will study how these two models violate the Wiedemann-Franz law, which states that for macroscopic ordinary metals the Lorenz ratio $L=K/GT$ is a constant equal to the Lorenz number $L_{0}=(k_{B}/e)^{2}(\pi^{2}/3)$ . Using Eqs. (48) and (58), we find that

[TABLE]

and

[TABLE]

where $L$ refers to the interacting case, and the expression of $K$ , for simplicity, has been computed at $\xi=0$ . In both cases the Wiedemann-Franz law is strongly violated: at $\Delta_{\text{min}}=0$ , the Lorenz ratio is exponentially smaller than $L_{0}$ thanks to $(\Delta E/k_{B}T)^{2}\exp{(-\Delta E/k_{B}T)}$ (this has been noticed in Ref. Tsaousidou2010, for the interacting model). In both cases the Lorenz ratio exponentially increases with $\Delta_{\text{min}}$ . In the non-interacting model the exponent is $2\Delta_{\text{min}}/(k_{B}T)$ (twice the interacting case), and the maximum value, achieved at $\Delta_{\text{min}}=\Delta E/2$ , is of the order of $L_{\text{NI}}\approx(\Delta E/k_{B}T)^{2}L_{0}$ . Interestingly, in the interacting case [Eq. (60)], when $|\Delta_{\text{min}}|>\Delta E$ , i.e. $N_{J}\neq 0$ , we find plateaus whose height increases with $N_{J}$ .

III.4 Three-terminal system

In this Section we consider the case of a three-terminal system, which allows to study the non-local transport coefficients and the influence of an additional terminal on the thermoelectric performance of the system. We focus on the simplest case when the couplings to the reservoirs are energy independent, that is, the rates $\Gamma_{\alpha}$ do not depend on $p$ , and we will consider an equidistant QD spectrum. Analytical expressions in the quantum limit for the power and the efficiency at maximum power can be obtained also in this case on the basis of the expressions written in Ref. Mazza2014, and of Eqs. (48) for the transport coefficients. All considerations made in this Section are valid for both interacting and non-interacting systems.

In a three-terminal setup with time reversal symmetry we have nine independent coefficients: three electrical conductances $G_{22}$ , $G_{33}$ , and $G_{23}$ , three thermopowers $S_{22}$ , $S_{33}$ , and $S_{23}$ , and three thermal conductances $K_{22}$ , $K_{33}$ , and $K_{23}$ . According to the expressions in Eqs. (48), valid in the quantum limit when the couplings to the leads are independent of energy, one finds that the local and non-local electrical (thermal) conductances are characterized by peaks (plateaus) located in the same positions as for the two-terminal case. The two local electrical (thermal) conductances $G_{22}$ , $G_{33}$ , ( $K_{22}$ , $K_{33}$ ), can have different heights if the coupling to reservoirs $2$ and $3$ are different ( $\Gamma_{2}\neq\Gamma_{3}$ ), while the non-local conductances, $G_{23}$ and $K_{23}$ , are negative. On the other hand, the local thermopowers $S_{22}$ and $S_{33}$ are equal and exhibit the same fine structure as in the two-terminal case. Moreover, the non-local thermopower $S_{23}$ vanishes, as a consequence of energy independent tunneling rates.

The power and efficiency of a three-terminal system are defined in Eqs. (19) and (20). For definiteness, let’s consider $T_{3}\geq T_{2}\geq T_{1}$ , so $\Delta T_{3}\geq\Delta T_{2}\geq 0$ . In general, Carnot’s efficiency cannot be written only in terms of the temperatures $T_{1}$ , $T_{2}$ , and $T_{3}$ , but it depends on the details of the system Mazza2014 . Nonetheless, if we fix the temperature of the hottest and coldest reservoir, that is, $T_{3}$ and $T_{1}$ , it can be shown that $\eta_{C}\leq\eta_{C}^{(2)}$ , where $\eta_{C}^{(2)}=1-T_{1}/T_{3}$ is the Carnot efficiency of the two-terminal system; the equal sign can be achieved when two reservoirs have the same temperatures. So adding a third terminal at an intermediate temperature cannot increase the maximum efficiency beyond the two terminal Carnot’s efficiency.

Also in the three terminal case the efficiency at maximum power cannot go beyond the linear response Curzon-Ahlborn efficiency $\eta_{CA}=\eta_{C}^{(2)}/2$ . Our aim is to maximize the power $P$ with respect to $\Delta\mu_{2}$ and $\Delta\mu_{3}$ , at given temperature differences. Then we will consider a fixed value of $\Delta T_{3}$ , and we will study the maximum power $P_{\text{max}}$ and the efficiency at maximum power $\eta(P_{\text{max}})$ varying $T_{2}$ between the fixed $T_{1}$ and $T_{3}$ . These calculations can be performed by writing the currents in terms of the temperature and electrochemical potential differences through the Onsager matrix $L_{ij}$ Mazza2014 ; in turn, $L_{ij}$ can be related to the transport coefficients. By rewriting Eqs. (43) as

[TABLE]

where

[TABLE]

and by defining

[TABLE]

we can write the currents as

[TABLE]

Note that the quantities $\mathcal{G}$ , $\mathcal{S}$ , $\mathcal{K}$ and $\mathcal{Z}T$ only depend on the properties of the QD, which can be interacting or non-interacting, and on the reference electrochemical potential $\mu$ ; they do not depend on the number of reservoirs nor on the tunneling rates.

The electrochemical potential differences that maximize $P$ at given reservoir temperatures can be written as

[TABLE]

Inserting these expressions into $P$ yields

[TABLE]

where $\mathcal{Q}=\mathcal{G}\mathcal{S}^{2}$ . The maximum power is thus an always positive quantity and, in the same manner as in the two-terminal case, see Eq. (21), it is proportional to $\mathcal{Q}$ and quadratic in the temperature differences. Furthermore, the properties of the QD and the chosen $\mu$ are all contained in the $\mathcal{Q}$ term, while the coupling to the reservoirs and the temperature differences are separately accounted for in the term between square parentheses in Eq. (68). The efficiency at maximum power instead is given by

[TABLE]

It is interesting to notice that in this equation the term before square parenthesis, which does not depend on $\Gamma_{\alpha}$ , is exactly equal to the two terminal efficiency at maximum power, see Eq. (22), since $\mathcal{Z}T=ZT$ for a two terminal system. The remaining part of Eq. (69) takes into account the particular temperatures and couplings to the three terminals. Furthermore the efficiency at maximum power, also in this 3 terminal system, only depends on $\mathcal{Z}T$ . There are three limiting cases we will first study: $\Delta T_{2}=0$ , $\Delta T_{2}=\Delta T_{3}$ and $\Delta T_{2}=\Delta T_{3}\Gamma_{3}/(\Gamma_{1}+\Gamma_{3})$ .

If $\Delta T_{2}=0$ , $T_{1}=T_{2}$ , so we have one hot reservoir at temperature $T_{3}$ and two cold ones at the same temperature. In this case we obtain

[TABLE]

Note that for a two-terminal system the maximum power, obtained by inserting Eq. (62) into Eq. (21), is given by

[TABLE]

Comparing Eq. (72) with Eq. (70), and Eq. (22) with Eq. (71), we notice that the 3 terminal system is formally equivalent to a 2-terminal system with temperature difference $\Delta T_{3}$ , with tunneling rate $\Gamma_{1}+\Gamma_{2}$ instead of $\Gamma_{1}$ , and $\Gamma_{3}$ instead of $\Gamma_{2}$ .

If $\Delta T_{2}=\Delta T_{3}$ , $T_{2}=T_{3}$ , so we have two hot reservoir at temperature $T_{3}$ and one cold reservoir at temperature $T_{1}$ . In this case we obtain

[TABLE]

As in the previous limiting case, this system behaves like a 2-terminal with temperature difference $\Delta T_{3}$ and with tunneling rate $\Gamma_{2}+\Gamma_{3}$ instead of $\Gamma_{2}$ .

If $\Delta T_{2}=\Delta T_{3}\Gamma_{3}/(\Gamma_{1}+\Gamma_{3})$ , reservoir 2 has an intermediate temperature such that $J^{h}_{2}=0$ . In fact this specific value of $\Delta T_{2}$ distinguishes the two regimes where $J^{h}_{2}>0$ and $J^{h}_{2}<0$ . In this case we obtain

[TABLE]

As in the other two limiting cases, this system behaves like a 2-terminal system, where reservoir 2 has been removed; this is to be expected because at this particular temperature, no heat flows through the second reservoir.

According to Eq. (69), all values of $\Delta T_{2}$ other than the three cases discussed above decrease the efficiency at maximum power with respect to the two-terminal case. The maximum power instead, at given tunneling rates, is increased with respect to the two-terminal case. In fact the maximum power is formally equal to that of a two-terminal system, coupled to the same QD, with increased tunneling rates. So if we have fixed values of $\Gamma_{1}$ , $\Gamma_{2}$ and $\Gamma_{3}$ , we achieve the largest maximum power by choosing $\Delta T_{2}=0$ if $\Gamma_{1}<\Gamma_{3}$ and $\Delta T_{2}=\Delta T_{3}$ if $\Gamma_{1}>\Gamma_{3}$ . In Fig. 5, $P_{\text{max}}$ and $\eta\left(P_{\text{max}}\right)$ from Eqs. (68) and (69) are plotted as a function of $\Delta T_{2}/\Delta T_{3}$ , choosing $\hbar\Gamma_{1}=0.02k_{B}T$ and $\hbar\Gamma_{2}=\hbar\Gamma_{3}=0.01k_{B}T$ . As we can see in panel (a), the power is maximum when $\Delta T_{2}=\Delta T_{3}$ (this result is expected since $\Gamma_{1}>\Gamma_{3}$ ). Furthermore the three maxima in panel (b) correspond to the three limiting cases previously studied, where $\eta\left(P_{\text{max}}\right)$ reaches the two-terminal performance.

IV Non-linear response

The linear response theory describes correctly the thermoelectric properties of bulk materials in most experimental conditions. However, as discussed for instance in Ref. bib:whitney, , non-linear effects are important in nanoscopic setups, since the temperature difference is applied across very small elements of the order of tens or hundreds of nanometers. As far as heat-to-work conversion is concerned, there is a practical reason to consider the non-linear response, namely efficiency and power output may increase with increasing temperature difference. Furthermore, for systems with time-reversal symmetry the efficiency at maximum power can overcome the limit of $\eta_{C}/2$ only beyond the linear response 444Note, however, that the Curzon-Ahlborn limit can be overcome also within linear response when time-reversal symmetry is broken, see Refs. Benenti2011, ; Balachandran2013, ; Brandner2013, ; Brandner2013b, ; Yamamoto2016, ..

In this Section we will consider a two-terminal QD system and discuss the numerical results obtained solving the kinetic equations (4) as discussed in Sec. II. We will focus our discussion on the thermoelectric properties, and on the efficiency and power produced by a QD-based heat engine. Let us define the charge current $J^{c}\equiv J^{c}_{2}=-J^{c}_{1}$ , thanks to charge current conservation, and the average reservoir temperature $\bar{T}=(T_{1}+T_{2})/2$ , which determines the typical thermal energy scale of the system beyond linear response (all energies will be given in units of $k_{B}\bar{T}$ ). Furthermore, $\Delta\mu\equiv\Delta\mu_{2}=eV$ , with $V$ applied voltage, $\Delta T\equiv\Delta T_{2}$ , and assume equidistant energy levels with spacing given by $\Delta E$ . In order to describe the potential drop between the QD and the two reservoirs, we will assume that the set of energy levels is shifted as $E_{p}(V)=E_{p}+(1-\theta_{0})eV$ , where $0\leq\theta_{0}\leq 1$ is the fraction of potential $V$ that drops over the tunnel barrier which couples reservoir 2 to the QD.

Regarding charge transport, we recall that in the linear response, within the quantum limit, the conductance exhibits peaks occurring at the dominant transition energies, i.e. when $\mu=\mu_{N}$ , of width given by $k_{B}\bar{T}$ . By applying a finite voltage bias $V$ , in the absence of a temperature difference ( $\Delta T=0$ ), the differential conductance, defined as

[TABLE]

exhibits the typical Coulomb diamond structure, with visible excited states, as a function of $V$ and $\mu$ , see for instance Ref. Vandersypen, .

IV.1 Non-linear Seebeck and Peltier coefficients

In the non-linear regime the thermopower (Seebeck coefficient) can be defined as follows

[TABLE]

i.e. as the ratio between the thermovoltage $V$ that develops as a result of a finite $\Delta T$ applied, at open circuit ( $J^{c}=0$ ).

In Fig. 6, $S$ is plotted as a function of $\mu$ for various values of $\Delta T/\bar{T}$ . The black solid curve (calculated for $\Delta T/\bar{T}=10^{-4}$ ) is the linear-response reference that is well approximated by the expression given in Eqs. (48). As discussed in Sec. III.1, the black solid curve presents main oscillations of period $\Delta E+2E_{C}$ , and a fine structure with a $\Delta E$ spacing bib:beenakker2 . Since we have chosen an equidistant energy spectrum, all curves share a number of features with the linear-response reference. Namely, i) $S$ crosses zero with positive slope at the main transition energies $\mu_{N}$ and is periodic with periodicity $\Delta E+2E_{C}$ (in Fig. 6 $\mu=290$ $k_{B}\bar{T}$ corresponds [see Eq. (49))] to $\mu_{N=3}$ ); ii) in the range of $\mu$ considered, $S$ is antisymmetric with respect to $\mu=290$ $k_{B}\bar{T}$ ; iii) $S$ vanishes for values of $\mu$ in the middle points between two dominant transitions $\mu_{N}$ and $\mu_{N+1}$ (in Fig. 6 such points are located at $\mu=230$ $k_{B}\bar{T}$ and $\mu=350$ $k_{B}\bar{T}$ ). Moreover, since we set $\theta_{0}=1/2$ , the linear increase of $S$ for $\mu\simeq\mu_{N}$ does not depend on the ratio $\Delta T/\bar{T}$ , i.e. it is well described by the linear-response proportionality coefficient $1/(|e|\bar{T})$ . Interestingly, such features (except for the fine structure oscillations) can be understood in terms of a non-interacting model (see App. D), which also explains the reduction of the negative slope of $S$ at the middle points as the ratio $\Delta T/\bar{T}$ increases.

Let us now discuss the behavior of $S$ when departing from the linear response regime. Fig. 6 shows that for all values of $\Delta T$ the thermopower deviates from the linear-response curve only for $\mu$ above 310 $k_{B}\bar{T}$ (or below 270 $k_{B}\bar{T}$ ). A sharp departure already occurs at $\Delta T/\bar{T}=0.5$ (red curve). This can be understood from the fact that, for $\mu>310$ $k_{B}\bar{T}$ , $S$ is of the order of $15~{}k_{B}/|e|$ which corresponds to a value of the thermovoltage ( $V=-7.5~{}k_{B}\bar{T}/|e|$ ) such that $|eV|\gg k_{B}T$ 555We have also checked that we are well beyond the linear response for the electric conductance.. Note that $\mu=310$ $k_{B}\bar{T}$ roughly corresponds to the first step of the fine structure in linear response. In particular, while the first step hardly moves by increasing $\Delta T$ from its position in linear response, the second step, occurring at $\mu=330$ $k_{B}\bar{T}$ in the linear response, shifts to a smaller value for increasing $\Delta T/\bar{T}$ , eventually disappearing or merging with the first step. This behavior may be attributed to the combination of the following two effects. On one hand, the thermovoltage $V$ , which determines the transport energy window, depends on $\mu$ and increases with $\Delta T$ according to the definition (78). On the other hand, an increase of $\Delta T/\bar{T}$ moves the lowest temperature ( $T_{1}$ ) towards absolute zero, thus sharpening the Fermi distribution function $f_{1}(E)$ . This last effect is also responsible for the sharpening of $S(\mu)$ as $\Delta T/\bar{T}$ increases. Note furthermore that the extremal values of $S$ decrease as $\Delta T/\bar{T}$ increases.

Let us now consider the non-linear Peltier coefficient defined, for a given voltage $V$ , as

[TABLE]

Our aim is to assess the failure of the Onsager reciprocity relation $\Pi=TS$ , which holds in the linear response regime. Beyond linear response, for a single-level non-interacting QD, one finds a “corrected” reciprocity relation, namely $\Pi+V/2=\bar{T}S$ in the case where $\theta_{0}=1/2$ (see App. D).

To single out the effect of interactions in a multi-level QD, in Fig. 7 we plot the ratio $r=(\Pi+V/2)/(\bar{T}S)$ as a function of $\mu$ for various values of $\Delta T/\bar{T}$ and $|e|V/(k_{B}\bar{T})$ ( $\Delta T/\bar{T}$ is used to compute $S$ , while $|e|V/k_{B}T$ to compute $\Pi$ ). Fig. 7 shows that the ratio $r$ departs significantly from 1, the linear response result (blue thin line), only far enough from the dominant transition energy $\mu_{N=3}=290$ $k_{B}\bar{T}$ . In particular, when $\Pi$ is in the linear response regime and $S$ is not (black solid curve) the strong deviations occurring for $\mu$ around the middle points between dominant transition energies, namely $\mu=230$ $k_{B}\bar{T}$ and $\mu=350$ $k_{B}\bar{T}$ , can be explained by a two-level non-interacting model (see App. D). However, the deviations occurring in the range of values of $\mu$ between 250 $k_{B}\bar{T}$ and 330 $k_{B}\bar{T}$ can be imputed to interaction effects. In the opposite case, where $S$ is in the linear response and $\Pi$ is not (red curve), the deviations of $r$ from 1 are entirely due to interaction effects and $r$ takes values between 0 and 1 in the entire range of values of $\mu$ . When both $S$ and $\Pi$ are beyond linear response the two behaviors discussed above coexist giving rise to the black dashed curve 666The black dashed curve is actually equal to the product of the solid black and the red curve. Indeed, due to the Onsager reciprocal relations in this product the linear response Peltier $\Pi$ (which appears in the ratio $r$ for the solid black curve) simplifies with the linear response $\bar{T}S$ (which appears in the ratio $r$ for the red curve)..

IV.2 Efficiency and output power

In this Section we consider the efficiency for heat-to-work conversion and output power in a two terminal system. Specifying Eq. (20) to a two terminal system where $\Delta T>0$ , we have that

[TABLE]

where $P>0$ is the output power, defined in Eq. (19), and $J^{h}_{2}>0$ is the heat current absorbed by the system. Apart from the system’s parameters, $\eta$ depends on $V$ and $\Delta T$ .

Let us first consider the maximum efficiency $\eta_{\text{max}}$ , obtained by maximizing the efficiency $\eta$ with respect to the applied voltage $V$ , at given $\Delta T$ . $\eta_{\text{max}}$ is plotted, normalized to $\eta_{C}$ , in Fig. 8(a) as a function of $\mu$ for different values of $\Delta T/\bar{T}$ . All plots show pairs of peaks close to $\mu_{N}$ , whose maximum is very close to $\eta_{C}$ , and secondary peaks of smaller height.

The solid black curve, relative to the linear response regime ( $\Delta T/\bar{T}=10^{-4}$ ), is related through Eq. (23) to the plot of $ZT$ [Fig. 3(b)]777Notice that in Fig 8 $\Delta E=20\,k_{B}T$ , while in Fig. 3 $\Delta E=10\,k_{B}T$ . discussed in Sec. III.2. For the black curve a pair of maxima approaching $\eta_{C}$ occur at $\mu=\mu_{N}\pm 2.40k_{B}\bar{T}$ , while $\eta_{\text{max}}$ vanishes at the dominant transition energies (and at the middle points between two dominant transition energies). Moreover, a fine structure of secondary peaks, with spacing $\Delta E$ , appears for intermediate values of $\mu$ . Moving away from the linear response, the main observation is that an increase of $\Delta T/\bar{T}>0.1$ produces only quantitative changes to the curves. As shown in Fig. 8(a), the main peaks of $\eta_{\text{max}}$ are still approximately located at $\mu=\mu_{N}\pm 2.40k_{B}\bar{T}$ and approaching the Carnot efficiency, while the peaks’ width reduces slightly with increasing $\Delta T$ , and the fine structure of the secondary peaks gets simply distorted.

Another important quantity in heat-to-work conversion is the maximum output power generated $P_{\text{max}}$ , which is obtained by maximizing the output power with respect to the applied voltage $V$ . It turns out that $P_{\text{max}}$ exhibits pairs of peaks approximately located at $\mu=\mu_{N}\pm 2.40k_{B}\bar{T}$ whose height increases approximately quadratically with $\Delta T$ , as long as $\Delta T/\bar{T}$ is not to close to 2. Interestingly, Fig. 8(b) shows that the maximum output power $P_{\text{max}}$ , when normalized to its peak value $P_{\text{peak}}$ , only very weakly depends on the ratio $\Delta T/\bar{T}$ . In particular, $P_{\text{max}}/P_{\text{peak}}$ is well approximated by the linear-response result, whose analytical expression is obtained by substituting Eq. (50) into Eq. (21).

The efficiency at maximum power $\eta(P_{\text{max}})$ can now be calculated by taking, for each value of $\mu$ , the value of $V$ which maximizes the power. $\eta(P_{\text{max}})$ is plotted, as a function of $\mu$ , in Fig. 8(c) for various values of the ratio $\Delta T/\bar{T}$ . By increasing such ratio starting from the linear response [solid black curve, related to the plot of $ZT$ in Fig. 3(b) through Eq. (22)] one finds that the peak values, again occurring approximately at $\mu=\mu_{N}\pm 2.40k_{B}\bar{T}$ , increase well above $\eta_{C}/2$ (the upper limit for the linear response). On the contrary, the efficiency at maximum power for values of $\mu$ away from $\mu_{N}$ , relative to the fine structure, decreases with increasing $\Delta T/\bar{T}$ beyond the linear response, but is only slightly different moving from $\Delta T/\bar{T}=1$ (red curve) to $\Delta T/\bar{T}=1.95$ (black dashed curve).

It is now interesting to compare the peak values of the efficiency at maximum power with various reference values, such as the Curzon-Alhborn (CA) efficiency bib:curzon $\eta_{\text{CA}}=1-\sqrt{1-\eta_{C}}$ and the Schmiedl-Seifert (SS) efficiency SchmiedlSeifert2008 $\eta_{\text{SS}}=\eta_{C}/(2-\eta_{C})$ 888This is actually the upper bound of the Schmiedl-Seifert efficiency, see Refs. SchmiedlSeifert2008, ; bib:esposito2, .. To do so we calculate the peak power, i.e. maximizing the power with respect to $V$ and $\mu$ , and plot the corresponding efficiency as a function of $\eta_{C}$ (determined by the temperature difference) in various situations, see Fig. 9. In particular, we consider the case of a QD with one doubly degenerate level with a finite charging energy (solid thick curves) and the case of a QD with two non-degenerate levels and zero charging energy (dashed curves). The parameters are chosen such that the two situations can be compared, namely $\Delta E=2E_{C}$ , i.e. the differential conductance consists of two peaks separated by the same electrochemical potential. Figure 9 (inset) shows the following general feature for small $\eta_{C}$ , i.e. in the linear response regime: $\eta(P_{\text{peak}})$ increases linearly with $\eta_{C}$ with slope determined by the value of $ZT$ [see Eq. (22)]. In particular, the two black curves (relative to $E_{C}=10$ $k_{B}\bar{T}$ and to $\Delta E=20$ $k_{B}\bar{T}$ ) virtually coincide, and are equal to the one for a single non-interacting level QD [see App. D and the first line of Eqs. (125)], since the parameters are such that $k_{B}\bar{T}\ll E_{C},\Delta E$ , where the transport is mostly accounted for by a single energy level. On the contrary, the two red curves (relative to $E_{C}=1.5$ $k_{B}\bar{T}$ and to $\Delta E=3$ $k_{B}\bar{T}$ ) differ by a large extent, with the interacting case (solid red curve) exhibiting larger efficiency at maximum power with respect to the associated non-interacting case (dashed red curve). Note that the efficiency at maximum power relative to the interacting case goes beyond the CA efficiency, when $\eta_{C}$ is larger than about 0.5, for all values of $E_{C}$ between 1.5 $k_{B}\bar{T}$ and 10 $k_{B}\bar{T}$ . Finally, we find that the SS efficiency is never overcome.

To complete the analysis, we show the correlation between efficiency and power corresponding to a few points (i.e. a few values of $\eta_{C}$ ) in the curve of Fig. 9. More precisely, Fig. 10 shows how the value of the power $P$ and the efficiency $\eta$ evolve by increasing the applied voltage $V$ from zero (where both $P$ and $\eta$ vanish) to the thermovoltage (where $P$ vanishes as a consequence of the fact that the charge current vanishes) 999For a given value of $\eta_{C}$ , we first find the values $\mu=\bar{\mu}$ and $V=\bar{V}$ leading to the peak value $P_{\text{peak}}$ for power. We then set $\mu=\bar{\mu}$ and vary the voltage $V$ from zero to the thermovoltage.. In particular, panel (a) refers to the case $\Delta E=20$ $k_{B}\bar{T}$ (dashed curves) and $E_{C}=10$ $k_{B}\bar{T}$ (solid curves), while panel (b) refers to the case $\Delta E=3$ $k_{B}\bar{T}$ (dashed curves) and $E_{C}=1.5$ $k_{B}\bar{T}$ (solid curves). We checked that in the linear response (when $\eta_{C}\ll 1$ ) the power reaches its maximum when the efficiency is nearly equal to $\eta_{C}/2$ . By increasing $\eta_{C}$ , for all the curves in the figure, both maximum power and efficiency at maximum power increase. For large values of $\eta_{C}$ (green curves), the efficiency remains close to the Carnot efficiency when $V$ is increased beyond the point of maximum power. The general feature is that in the interacting case (solid curves) the power is much larger than in the associated non-interacting situation (dashed curves). When $\Delta E$ and $E_{C}$ are of the same order as $k_{B}\bar{T}$ [Fig. 10(b)], both the maximum power and the efficiency at maximum power are increased in the interacting case as compared with the non-interacting case. A remarkable property of both regimes discussed in panels (a) and (b) is that in the strongly nonlinear regime, the maximum power is obtained for values where the efficiency is high and close to the maximum efficiency.

We finally discuss a peculiarity of the maximum power output considering a doubly degenerate energy level in the interacting case 101010The effect of degenerate energy levels on the electron and thermal conductance and on thermopower of a QD was considered in Ref. bib:zianni_2008, . (as we will point out, some considerations are also valid in the non degenerate case).

Fig. 11 shows the maximum power, normalized with respect to the peak value, as a function of $\mu$ for different values of $\eta_{C}$ including the linear-response case (parameter values are the same as for Fig. 9, with $E_{C}=10$ $k_{B}\bar{T}$ ). The first remarkable feature is that the two peaks around a dominant transition energy have different heights (here $\mu_{N=1}=10$ $k_{B}\bar{T}$ and $\mu_{N=2}=30$ $k_{B}\bar{T}$ ); in the absence of degeneracy, within the quantum limit all the peaks have the same height (see Fig. 3). More precisely, see Eq. (81), the external peaks, displaced of $2.53k_{B}\bar{T}$ from $\mu_{N}$ are higher with respect to the internal peaks, displaced of $2.32k_{B}\bar{T}$ from $\mu_{N}$ , whose height is almost equal to the non-degenerate (interacting) case. The second feature is that all curves, apart from the case of very large $\eta_{C}$ (blue curve), are well approximated by the linear-response quantum limit expression (black dashed curve in Fig. 11)

[TABLE]

The plus (minus) sign in Eq. (81) is to be taken when $N_{\text{min}}$ is even (odd). Eq. (81) allows us to find that level degeneracy gives rise to an enhancement of output power of about 1.77 times respect to the non-degenerate case, independently of the parameter values (see also Ref. kuo2016, ). Eq. (81) also makes clear that the origin of the asymmetry for $\Delta_{\text{min}}\to-\Delta_{\text{min}}$ and of the difference in peaks’ height is the term $\pm\log 2/2$ occurring in the presence of degeneracy. The case of $E_{C}=1.5$ $k_{B}\bar{T}$ (not shown), which is not within the quantum limit, produces a much more asymmetric behavior. We have verified that also in the non degenerate case the analytic formula derived in the linear response regime [see Eq. (50)] well describes the maximum power also beyond the linear response regime.

V Conclusions

We have studied the thermoelectric properties of a multi-level interacting QD in the sequential tunneling regime, in a multi-terminal setup, both in the linear response regime and beyond. In particular, we have

•

generalized the sequential tunneling method put forward by Beenakker in Refs. bib:beenakker1, and bib:beenakker2, to a multi-terminal configuration and set the range of validity of the expressions for the charge and heat currents in the linear response regime;

•

found simple analytic formulas for the multi-terminal transport coefficients in the low temperature limit;

•

found simple analytic formulas for the power factor $Q$ and the figure of merit $ZT$ in the low temperature limit for a two-terminal setup;

•

found that $Q$ and $ZT$ can be simultaneously maximized for suitable values of the electrochemical potential;

•

found that Coulomb interactions can dramatically enhance $ZT$ by suppressing the thermal conductance;

•

found that both the interacting an non-interacting models strongly violate the Wiedemann-Franz law;

•

found analytic expressions for the maximum power and for the efficiency at maximum power in a three terminal setup;

•

investigated the nonlinear Seebeck and Peltier coefficients in a two-terminal setup, identifying features of the breakdown on the Onsager reciprocity relation;

•

computed numerically the maximum efficiency, the maximum power, and the efficiency at maximum power in the non-linear regime, finding optimal system parameters for heat-to-work conversion such that the efficiency at maximum power can go beyond Curzon-Alhborn’s efficiency;

•

compared the case of a doubly degenerate level with interaction and the case of two non-degenerate levels without interaction finding that the interacting case enhances the power output and, especially when charging energy and level spacing are of the order of the thermal energy, it increases the efficiency at maximum power that can go beyond Curzon-Alhborn’s efficiency;

•

found that the non-linear maximum power is well approximated by the analytic linear response expression;

•

found that QDs with degenerate energy levels and Coulomb interactions achieve higher efficiency and output power than non-degenerate QDs; in particular the maximum power is enhanced almost of a factor 2;

•

calculated the transport coefficients for a non-linear, non-interacting QD with 1 and 2 energy levels (App. D);

•

found approximate analytic expressions for the maximum power and efficiency at maximum power for a non-linear, non-interacting QD with 1 energy level (App. D).

The multi-terminal formalism developed in this paper and the expressions we have obtained for charge and heat currents, transport coefficients, power and efficiency could be used to design and analyze experimental data. Extensions of the studies presented in this paper could include level spacings different from the equidistant and regimes beyond the quantum limit. Finally, a comprehensive description of the thermoelectric properties and performance of a QD should assess the role of quantum coherence going beyond the sequential tunneling limit and the relevance of phonon contribution to heat transport.

Acknowledgements

We would like to acknowledge fruitful discussions with Stefano Roddaro. This work has been supported by the SNS internal projects “Thermoelectricity in nanodevices” and “Non-equilibrium dynamics of one-dimensional quantum systems: From synchronization to many-body localization”, by the EU project ThermiQ and by the COST Action MP1209 “Thermodynamics in the quantum regime”. G.B. acknowledges support by INFN through the project “QUANTUM”.

Appendix A The kinetic equations always allows a non-trivial solution.

To prove this statements, let us consider the kinetic equations as written in Eq. (11). This equation is made up of a sum over $p$ of the following expression:

[TABLE]

which is a function of a generic configuration $\{n_{i}\}$ . Let us sum this expression for the two particular configurations: $(\{n_{i}\},n_{p}=1)$ and $(\{n_{i}\},n_{p}=0)$ . When $n_{p}=1$ , the first term in round parenthesis gives a plus sign, while when $n_{p}=0$ , it gives a minus sign. The rest of the expression does not depend on $n_{p}$ , since $\widetilde{N}=\sum\nolimits_{i\neq p}n_{i}$ , thus the sum over the above configurations exactly vanishes. Now let us go back to Eq. (11), and let us sum over the two configurations $(\{n_{i}\},n_{k}=1)$ and $(\{n_{i}\},n_{k}=0)$ , $k$ being a given index. According to the argument given above, the term in the sum where $p=k$ vanishes, and we obtain

[TABLE]

Thus by summing over a given occupation number $n_{k}=0,1$ we have removed the case $p=k$ in the sum over $p$ . If we now sum over all occupation numbers, we will remove all terms from the sum, yielding zero 111111By summing over all occupation numbers except for a fixed $n_{p}$ , we could demonstrate that the LBEs derive from the kinetic equations.. The sum over all occupation numbers $\{n_{i}\}$ is the sum of all $2^{L}$ equations of the kinetic equations: this demonstrates that any equation in the kinetic equations is linearly dependent from the other ones. Thus the matrix $M$ defined by the kinetic equations has a null space of dimension at least 1, since we demonstrated that the rows of $M$ are not linearly independent. Furthermore, if we perform the same sum over all occupation numbers to the time dependent kinetic equations, given in Eq. (4), we find that:

[TABLE]

This is an obvious but important property that says that the probability normalization does not depend on time.

Appendix B The DBEs are not consistent in general

In general, if $E_{C}\neq 0$ , the DBEs are not consistent. This means that no set of $P(\{n_{i}\})$ exists that can simultaneously satisfy all the DBEs. In App. (C) we will discuss which conditions guarantee their consistency within the linear response regime.

Here let us demonstrate their inconsistency in the special case $L=2$ . Let us start from Eq. (15) and consider non null temperature and electrochemical potential differences. We will show that these equations form an over-complete set for $P(\{n_{i}\})$ that in general does not allow any non null solution. Since $L=2$ , we have $2^{L}=4$ unknown probabilities, and the number of DBEs is $2^{L-1}2=4$ . In this case, the DBEs can be represented in matrix form as follows

[TABLE]

where $P_{n_{2}n_{1}}\equiv P(n_{2},n_{1})$ and the coefficients are defined in Eqs. (12) and (13). In order to show that this linear algebra problem does not allow a non null solution, we will show that the determinant $\det(M_{D})$ is in general not zero. We have

[TABLE]

where $\Gamma_{\text{tot}}(p)=\sum_{\alpha}\Gamma_{\alpha}(p)$ . It is pretty clear that since $\Gamma_{\alpha}(p)$ , $E_{p}$ and $E_{C}$ are arbitrary, this determinant cannot be in general zero. For instance, choosing a two terminal system with $\Gamma_{1}(p)=\Gamma_{2}(p)=\Gamma$ , $\Delta\mu_{2}=4k_{B}T$ , $E_{1}=1/2k_{B}T$ , $E_{2}=3/2k_{B}T$ , $\mu=\Delta T_{2}=0$ we obtain $\text{det}(M_{D})/\Gamma^{4}\simeq 0.20$ . It is interesting to notice that at equilibrium $(\Delta\mu_{2}=\Delta T_{2}=0)$ , the DBEs are all exactly satisfied by the gran canonical distribution. Furthermore, in the non-interacting limit $E_{C}=0$ , the coefficients $A_{\widetilde{N},p}$ and $B_{\widetilde{N},p}$ depend only on $p$ , so that we can drop the $\widetilde{N}$ argument. Thus $\text{det}\left(M_{D}\right)=A_{2}A_{1}B_{1}B_{2}-A_{1}A_{2}B_{2}B_{1}\equiv 0$ .

The proof can be extended to any number $L$ of levels as follows. We rewrite the DBEs, Eq. (15), as

[TABLE]

This equation has the same form as the LDBEs (89) which we will consider in App. C, where the unknown probabilities are replaced by the logarithm of the probabilities, and where

[TABLE]

As we shall show in App. C, these equations are consistent if $\delta_{p}(\widetilde{N})$ satisfies property (91). This is the case if $E_{C}=0$ , since $A_{\widetilde{N},p}$ and $B_{\widetilde{N},p}$ are then independent of $\widetilde{N}$ . Then condition (91) is trivially fulfilled with $c=0$ and the DBEs are consistent.

To summarize, as we will see in Appendix C, the DBEs are not consistent in general, but they are if we set $E_{C}=0$ or $\Delta\mu_{\alpha}=\Delta T_{\alpha}=0$ for all $\alpha$ . Furthermore, within the linear response regime, they are valid in many more cases, for example if the tunneling rates are proportional, i.e. $\Gamma_{\alpha}(p)=k_{\alpha}\Gamma_{1}(p)$ for $\alpha=2,3,\dots,{\cal N}$ .

Appendix C Conditions of validity of the linearized DBEs

In this Appendix we will assess under which conditions the linearized DBEs are consistent. We start by rewriting Eq. (32) in the following way:

[TABLE]

where

[TABLE]

We will prove that the linearized DBEs, written in Eq. (89), are consistent if $\delta_{p}(\widetilde{N})$ satisfies the property

[TABLE]

where $c$ is a constant that does not depend on $p$ , $N$ or $M$ . This statement will be explicitly proven for a two energy level system, then it will be extended to $L$ energy levels by induction.

Property (91) is in general satisfied if the tunneling rates are proportional, i.e. $\Gamma_{\alpha}(p)=k_{\alpha}\Gamma_{1}(p)$ , or if $E_{C}=0$ , or if $\Delta T_{\alpha}=0$ for all $\alpha$ . Furthermore, in a three terminal system, property (91) is satisfied also if $\Delta T_{2}=0$ and $\Gamma_{3}(p)=k\Gamma_{\text{tot}}(p)$ , or if $\Delta T_{3}=0$ and $\Gamma_{2}(p)=k\Gamma_{\text{tot}}(p)$ . These conditions can be generalized to a generic ${\cal N}$ -terminal system by requiring that

[TABLE]

is independent of $p$ .

C.1 Two energy level system ( $L=2$ )

We will prove that if property (91) is satisfied, then the linearized DBEs allow a solution. Eq. (89) represents a linear algebra problem for $\psi(\{n_{i}\})$ . For $L=2$ , $\psi(\{n_{i}\})=\psi(n_{2},n_{1})\equiv\psi_{n_{2},n_{1}}$ . Let $\vec{\psi}$ be the vector $(\psi_{00},\psi_{01},\psi_{10},\psi_{11})$ and let $B$ be the corresponding matrix such that Eq. (89) can be written as

[TABLE]

Matrix $B$ has a null space of dimension 1, generated by $\vec{\psi}=(1,1,1,1)$ , thus it is not invertible. This vector actually represents the equilibrium distribution: when $\Delta T_{\alpha}=\Delta\mu_{\alpha}=0$ , $\delta_{p}(\widetilde{N})=0$ , so $\vec{\psi}=(1,1,1,1)$ satisfies Eq. (93). Since $P=P_{eq}(1+\psi)$ , $P=2P_{eq}$ . Normalizing the probabilities yields $P\equiv P_{eq}$ , so we have demonstrated that the equilibrium distribution is in fact given by the grand canonical distribution. Eq. (93) will allow a solution if and only if vector $\vec{\delta}$ belongs to the image of matrix $B$ . The dimension of the image of matrix $B$ is 3, so there is a one dimensional space orthogonal to the image of $B$ that cannot be obtained by linear combinations of $B$ ’s columns. Vector

[TABLE]

is orthogonal to the columns of $B$ , so it is a generator of this one dimensional space. Thus a solution $\vec{\psi}$ will exist if and only if vector $\vec{\delta}$ only belongs to the image of $B$ , thus it cannot have a projection on $\vec{v_{0}}$ . The projection is zero when

[TABLE]

which is satisfied using property (91).

As we have shown in Sec. III, Eq. (28) for the linearized charge current is correct in general, instead the linearized heat current given in Eq. (33) is correct only if the linearized DBE are consistent. We verified numerically all these statements computing exactly, for small temperature and electrochemical potential differences, the charge and heat currents using both the DBEs and the kinetic equations. As demonstrated by Fig. 12, which shows the heat current as a function of $\mu$ for a choice of tunneling rates that are not proportional to each other, a particular subset of DBEs leads to an incorrect result (red dashed curve) which differs from the one obtained using the kinetic equations (black solid curve). On the other hard, there is no difference when plotting the charge current using the same parameters as in Fig. 12.

C.2 $L$ energy level system

We extend the previous demonstration to a generic system with $L$ energy levels. We will use the following notation and conventions:

•

$\{n_{i}\}$ indicates a generic set of occupation numbers, conventionally ordered the following way: $\{n_{L},n_{L-1},\dots,n_{2},n_{1}\}$ .

•

$(n_{p}=1,\{n_{i}\})$ is a shorthand notation for: $\{n_{L},\dots,n_{p+1},1,n_{p-1},\dots,n_{1}\}$ .

•

$(n_{p}=0,\{n_{i}\})$ is a shorthand notation for: $\{n_{L},\dots,n_{p+1},0,n_{p-1},\dots,n_{1}\}$ .

•

$(n_{p}=1,\{0\})$ is a shorthand notation for: $\{0,\dots,0,1,0,\dots,0\}$ , where the $1$ is relative to $n_{p}=1$ .

•

$(n_{p}=0,\{0\})$ is equivalent to: $\{0,\dots,0\}$ .

•

$\vec{\psi}_{i}$ is the component of vector $\vec{\psi}$ with index $i$ .

•

Each linearized DBE in Eq. (89) can be uniquely defined by specifying $p$ and $\{n_{i}\}$ , so $[p|\{n_{i}\}]$ is a shorthand notation to indicate that given equation.

•

Let $A$ be a matrix written in block form, for example:

[TABLE]

where $B$ , $C$ , $D$ and $E$ are matrices. By “block row” we intend the set of rows in the block. For example, the first block row of $A$ is $(B,0)$ , the second block row of $A$ is $(0,C)$ , and so on.

In order to represent the linearized DBE, given in Eq. (89), in matrix form, we will put the $2^{L}$ unknowns $\psi(\{n_{i}\})$ into a vector as follows: $\vec{\psi}_{I(\{n_{i}\})}\equiv\psi(\{n_{i}\})$ , where $I(\{n_{i}\})=\sum\limits_{i=1}^{L}n_{i}2^{i-1}$ .

Let us define $B(L)$ as the $L2^{L-1}\times 2^{L}$ matrix such that Eq. (89) can be represented as follows:

[TABLE]

where $\vec{\delta}(L)$ is a vector containing the various values of $\delta_{p}(\widetilde{N})$ . Each row of $B(L)$ contains zeros except for a $+1$ and a $-1$ in the position corresponding to $\psi(n_{p}=0,\{n_{i}\})$ and $\psi(n_{p}=1,\{n_{i}\})$ . Let’s notice that $B(L)$ has the following property:

[TABLE]

This property is trivial since, as we just discussed, each row of $B(L)$ only has a $+1$ and a $-1$ .

C.2.1 Recursive decomposition

We want to write $B(L)$ in terms of $B(L-1)$ . We will thus order the rows in $B(L)$ , and accordingly the elements of $\vec{\delta}(L)$ , starting from those LDBE $[p|\{n_{i}\}]$ where $p<L$ , then we will add the ones involving level $L$ . In particular, we will first put equations $[p|(n_{L}=0,\{n_{i}\})]$ with $p<L$ , then we will add equations $[p|(n_{L}=1,\{n_{i}\})]$ with $p<L$ , and at last we will add equations $[L|\{n_{i}\}]$ , which are the ones involving level $L$ . For example, $B(2)$ will be written as in Eq. (93).

This ordering allows us to relate $B(L)$ to $B(L-1)$ . The first set of equations $[p|(n_{L}=0,\{n_{i}\})]$ , with $p<L$ , represents all possible LDBE between elements of $\vec{\psi}$ with $n_{L}=0$ , so they are equivalent to the LDBE for $L-1$ levels. Since $n_{L}=0$ , according to our indexing convention, these rows involve all indexes of $\vec{\psi}_{i}$ with $0\leq i<2^{L-1}$ , which corresponds to the first half of vector $\vec{\psi}$ . Instead, the equations $[p|(n_{L}=1,\{n_{i}\})]$ , with $p<L$ , represents all possible LDBE between elements of $\vec{\psi}$ with $n_{L}=1$ , so also these equations are equivalent to the LDBE for $L-1$ levels. Since $n_{L}=1$ , according to our indexing convention, these rows involve all indexes of $\vec{\psi}_{i}$ with $i\geq 2^{L-1}$ , which corresponds to the second half of vector $\vec{\psi}$ . At last, equations $[L,\{n_{i}\}]$ relate components of $\psi$ where only the occupation number $n_{L}$ is changed. Since $I(n_{L}=1,\{n_{i}\})=I(n_{L}=0,\{n_{i}\})+2^{L-1}$ , these equations relate indexes of $\vec{\psi}$ that are distant $2^{L-1}$ , thus they mix elements between the first and second half of $\vec{\psi}$ . We will thus have the following block representation:

[TABLE]

where $I_{d}(L)$ is the $2^{L}\times 2^{L}$ identity matrix. Since we trivially have that $B(1)=(1,-1)$ , Eq. (99) can be used to define $B(L)$ recursively, yielding a precise row ordering. For example, as we can see in Eq. (93), $B(2)$ can be obtained by applying Eq. (99) to $B(1)$ .

Also $\vec{\delta}(L)$ allows a decomposition in terms of $\vec{\delta}(L-1)$ . The first block row in Eq. (99), $(B(L-1),0)$ , corresponds to equations $[p|(n_{L}=0,\{n_{i}\})]$ with $p<L$ . Recalling Eq. (90), we will notice that $\delta_{p}(\widetilde{N})$ depends on level $L$ only through:

[TABLE]

since $n_{L}=0$ . So in the first block row we have that $\vec{\delta}(L)=\vec{\delta}(L-1)$ . The second block row, $(0,B(L-1))$ , corresponds to equations $[p|(n_{L}=1,\{n_{i}\})]$ with $p<L$ . So this time the presence of level $L$ will change the values of $\widetilde{N}$ the following way:

[TABLE]

since $n_{L}=1$ . So in the second block row $\vec{\delta}(L)$ is given by $\vec{\delta}(L-1)$ replacing $\widetilde{N}$ with $\widetilde{N}+1$ . Using property (91), we have that:

[TABLE]

where $c$ is a constant that does not depend on $p$ or $\widetilde{N}$ . So in the second block row we will have that $\vec{\delta}(L)=\vec{\delta}(L-1)+c$ . The last row $(I_{d}(L-1),-I_{d}(L-1))$ in Eq. (99) corresponds to equations $[L|\{n_{i}\}]$ , which involve level $L$ . So let us define $\vec{\partial}(L)$ as the vector that corresponds to $\vec{\delta}(L)$ in the third block row.

Using these observations and Eq. (99), we can rewrite Eq. (97) using the following decomposition:

[TABLE]

This decomposition will be the main tool to perform a demonstration by induction.

C.2.2 The independent part of $B(L)$

Among the $L2^{L-1}$ rows of $B(L)$ , we will explicitly show how to extract $2^{L}-1$ linearly independent rows; we will later show that all other rows can be obtained by linear combinations. Let’s define $\widetilde{B}(L)$ as the matrix containing only the rows of $B(L)$ corresponding to equations $[p|\{n_{i}\}]$ with $n_{i}=0$ for $i<p$ ; these are $2^{L}-1$ rows. We will now prove by induction that all the rows of $\widetilde{B}(L)$ are linearly independent.

Since $\widetilde{B}(L)$ is made by selecting rows from $B(L)$ , we can decompose it in terms of $\widetilde{B}(L-1)$ just like we did for $B(L)$ in Eq. (99). An interesting characteristic of $\widetilde{B}(L)$ is that it only has one equation in the block $[L|\{n_{i}\}]$ involving $\psi(0,0,\dots,0)$ and $\psi(1,0,\dots,0)$ . So the decomposition becomes:

[TABLE]

We can use Eq. (104) to define $\widetilde{B}(L)$ recursively noticing that $\widetilde{B}(1)=B(1)=(1,-1)$ .

Now we can prove by induction over $L$ the following statement: all the rows of $\widetilde{B}(L)$ are linearly independent. This statement is obvious for $L=1$ and $L=2$ . We will now assume that our statement is valid for $L-1$ , and we will show that this implies that it is valid for $L$ . The first block row in Eq. (104) contains rows that are linearly independent by hypothesis. The same is true for the second block row. Also putting together the first and second block row yields $2^{L}-2$ independent rows, since the space generated by the respective linear combinations are clearly orthogonal thanks to the block decomposition. So we have to show that the last row cannot be written as a linear combination of previous rows. We will make a proof by contradiction. Let’s assume that the last row $\vec{v}=(1,0,\dots,0|-1,0,\dots,0)$ can be written as a linear combination of a vector belonging to the space generated the first block row, $\vec{w_{0}}$ , and of one belonging to the second block row, $\vec{w_{1}}$ :

[TABLE]

If we define $\vec{n}=(1,\dots,1|0,\dots,0)$ , we will have that:

[TABLE]

The first equation is zero because, as we have shown in Eq. (98), $B(L)$ has the property that $B(L)(1,\dots,1)^{T}=0$ , so it is true also for $\widetilde{B}(L)$ . The second equation is zero because $\vec{n}$ has zeroes on all non null components of $\vec{w_{1}}$ . So multiplying each sides of Eq. (105) by $\vec{n}$ yields $\vec{n}\cdot\vec{v}=\vec{n}\cdot(\vec{w_{0}}+\vec{w_{1}})=0$ , which is absurd because $\vec{n}\cdot\vec{v}=1$ . Thus we have completed the proof by induction.

C.2.3 Consistency of the linearized DBEs

We will now show that if property (91) is valid, the linearized DBEs are consistent and allow a solution for $\psi$ given by a one dimensional space; the unique solution can then be found by imposing the normalization condition on the probabilities. We will thus prove by induction over $L$ the following statement: if $\delta_{p}(N)$ satisfies property (91), all equation in the linear algebra system in Eq. (103) can be written as linear combinations of the equations involving the rows of $\widetilde{B}(L)$ . Since there are $2^{L}-1$ independent rows in $\widetilde{B}(L)$ , by virtue of Rouché-Capelli’s theorem we will have a solution given by a $2^{L}-(2^{L}-1)=1$ dimensional space. This statement is trivially true for $L=1$ , since $B(1)=\widetilde{B}(1)$ , and we have proved it explicitly for $L=2$ , so we will show that if it’s valid for $L-1$ , then it’s valid for $L$ .

Let’s consider the subset of rows of Eq. (103) where $B(L)=\widetilde{B}(L)$ . Among these equations, we will denote by “set 0” those belonging to the first block row, and “set 1” those belonging to the second block row. As we can see from Eq. (104), this leaves out one equation in the last block row, which is:

[TABLE]

By induction hypothesis, all equations in the first block row of Eq. (103) can be written as linear combinations of set 0. Furthermore, if $\delta_{p}(N)$ satisfies property (91), also $\delta_{p}(N)+c$ satisfies the same property, so by induction hypothesis, also all equations in the second block row of Eq. (103) can be written as linear combinations of set 1. So any row in the first and second block row in Eq. (103) can be written as linear combinations of those equations involving $\widetilde{B}(L)$ .

We are left to show that equations in the last block row of Eq. (103) can be written as linear combinations of equations involving $\widetilde{B}(L)$ . Thanks to what we have just demonstrated, this is is equivalent to showing that the equations in the last block row are linear combinations of the first two block rows of Eq. (103), and of Eq. (107). The equations of the last block row of Eq. (103) are of the form

[TABLE]

Given any fixed set of $\{n_{i}\}$ , which fixes $\widetilde{N}=\sum_{i\neq L}n_{i}$ , we will now show how to obtain Eq. (108) by relating $\psi(n_{L}=0,\{n_{i}\})$ to $\psi(n_{L}=1,\{n_{i}\})$ only using the equations of the first two block rows of Eq. (103) and Eq. (107). We can first relate $\psi(n_{L}=0,\{n_{i}\})$ to $\psi(n_{L}=0,\{0\})$ using the equations in the first block row $\widetilde{N}$ times, each time changing a non null $n_{i}$ to zero. We will obtain the following result:

[TABLE]

where $p_{i}$ is the energy level whose occupation number has been changed to zero at step $i$ ; the argument $\widetilde{N}-i$ is given by the fact that at each step we use a linearized DBE where $\sum_{i\neq p_{i}}n_{i}$ decreases by one. Using Eq. (107) we obtain

[TABLE]

Now we can relate $\psi(n_{L}=1,\{0\})$ to $\psi(n_{L}=1,\{n_{i}\})$ using the equations of the second block row to set each non null $n_{i}$ back to $1$ . If we do this in the exact same order as we did in the previous step, we will obtain

[TABLE]

The “new” term has the same set of $\{p_{i}\}$ , and the argument $i$ is given by the fact that $\sum_{i\neq p_{i}}n_{i}$ increases by one at each step, and it starts from $1$ since now $n_{L}=1$ . Eq. (108) is of the same form, therefore the equations on the last block row can be expressed as linear combinations of $\widetilde{B}(L)$ if and only if they represent the same equation, which requires:

[TABLE]

Using property (91), this last equation reduces to $\widetilde{N}=\widetilde{N}$ , concluding the proof by induction.

Appendix D Non-interacting QD: Non-linear effects

In this Appendix we will study the transport coefficients, maximum and peak efficiency and efficiency at peak power of a two-terminal non-interacting QD. We will first consider a single-level QD and then a two-levels QD (with energy-independent tunneling rates).

Let us first consider a QD with a single energy level $E_{1}$ . In this case, the kinetic equations with $E_{C}=0$ reduces to a single simple equation, so it is possible to compute the currents. Using the two-terminal notation introduced in Secs. III.2 and IV we obtain

[TABLE]

where

[TABLE]

with the Fermi function $f(x)=[1+\exp(x)]^{-1}$ . Note that the same result could be obtained using the Landauer-Büttiker scattering formalism with a narrow in energy single-level QD transmission probability. Using Eq. (113), we can compute

[TABLE]

Note that the value of $\theta_{0}$ can be determined by measuring $G$ .

We can also compute $S$ using Eq. (113). The condition $J^{c}=0$ is satisfied when the arguments of the two Fermi distributions $f_{1}$ and $f_{2}$ are equal, so we obtain

[TABLE]

where

[TABLE]

is the average reservoir temperature, weighed with $\theta_{0}$ . As we can see, the slope of $S(\mu)$ beyond the linear response regime is strongly determined by $\theta_{0}$ . In the linear response regime $T^{*}\simeq T$ , where $T$ is the average reservoirs’ temperature, so once again we obtain the result we obtained in the quantum limit linear response regime of an interacting QD neglecting the fine structure oscillations, due to other energy levels.

It is now interesting to study the influence of a second energy level on $S$ . In Fig. 13 we show the comparison between the thermopower of a non-interacting two energy level QD (dashed lines) and the multilevel interacting QD considered in Fig. 6 (solid lines), for different values of $\Delta T$ . For the sake of comparison, we choose the value for the two energy levels $E_{1}$ and $E_{2}$ to match the two dominant transition energies of the multilevel interacting QD, in particular $E_{1}=\mu_{N=3}=290$ $k_{B}\bar{T}$ and $E_{2}=\mu_{N=4}=410$ $k_{B}\bar{T}$ .

On one hand, the non-interacting thermopower shows a linear dependence on $\mu$ when $\mu$ is close to $E_{1}$ or $E_{2}$ , whose slope is independent of $\Delta T$ (only because $\theta_{0}=1/2$ ) and is the same as the interacting QD when $|\Delta_{\text{min}}|<\Delta E$ . On the other, when $\mu$ is close to the middle point $\mu^{*}=(E_{1}+E_{2})/2$ between $E_{1}$ and $E_{2}$ , the thermopower depends linearly on $\mu$ with a negative slope ( $S\simeq\alpha\delta\mu$ , with $\delta\mu$ being a small displacement with respect to $\mu^{*}$ ) which depends on the temperatures as

[TABLE]

where $\Delta\varepsilon=E_{2}-E_{1}$ . Eq. (118) has been found by expanding the charge current around $\mu=\mu^{*}$ . Once again, the slope of the interacting and non-interacting QDs are the same in this region. In between instead the fine structure oscillations of the interacting QD create a substantial difference that causes the maximum value of $S$ to decrease.

At last we can compute the Peltier coefficient for a single energy level QD using Eq. (113):

[TABLE]

Surprisingly, we would have obtained the same result also without setting $\Delta T=0$ . On the other hand, it strongly depends on $V$ : increasing the voltage shifts $\Pi$ by $\theta_{0}V$ , so measuring this shift would allow to measure $\theta_{0}$ . Furthermore, if we compute $\Pi$ in the linear response regime, setting $V=0$ in Eq. (119), we can explicitly verify that $\Pi=TS$ is respected. If we set $\theta_{0}=0$ , we see that $\Pi=T^{*}S$ , even beyond the linear response regime, for arbitrary temperature and voltage differences. Interestingly, in the presence of two levels the Peltier coefficient beyond linear response $\Pi$ still satisfies the relation $(\Pi+\theta_{0}V)=TS_{\text{lin}}$ , where $S_{\text{lin}}$ is the linear-response thermopower (black dashed curve in Fig. 13), so that $(\Pi+\theta_{0}V)$ does not depend on $V$ .

The linearization of $G$ and $S$ for a single-level non-interacting QD yields the same results obtained for a multilevel interacting QD, in the quantum limit linear response regime, if we restrict $|\Delta_{\text{min}}|<\Delta E$ . We could expect this result since considering a single energy level is intuitively equivalent to sending $\Delta E$ and $E_{C}$ to infinity in the interacting QD. On the contrary, a single level QD model cannot be used to estimate $K$ since it predicts $K=0$ [in fact, from Eq. (113) we have that $J^{h}_{2}\propto J^{c}$ , so $K$ , caluclated at $J^{c}=0$ , vanishes].

We now want to study the power and efficiency of a single-level QD. Considering $\Delta T=T_{2}-T_{1}>0$ and inserting Eq. (113) into the definitions in Eqs. (19) and (20), we can write the power and efficiency as

[TABLE]

We will consider a fixed temperature difference $\Delta T=T_{2}-T_{1}>0$ , and a variable $V$ such that the system behaves as a heat engine ( $P>0$ ). The power is positive when $V\in[0,V_{\text{stop}}]$ , where $V_{\text{stop}}$ is the non zero voltage that creates a null charge current. Imposing this condition we find

[TABLE]

Without loss of generality, we can specify our analysis to the region where $\mu<E_{1}$ , $E_{1}$ being the energy of the single energy level; thus $\bar{\Delta}_{\text{min}}>0$ , so $V_{\text{stop}}$ is positive, and our system will behave as a heat engine when

[TABLE]

Let us now discuss the peak efficiency of the system. We have to maximize Eq. (121) with respect to $\bar{\Delta}_{\text{min}}$ and $V$ , at fixed $\Delta T>0$ , respecting $\bar{\Delta}_{\text{min}}>0$ and Eq. (123). Since $\eta$ is a growing function of $V$ for $\bar{\Delta}_{\text{min}}>0$ , $\eta$ will be maximum when computed at the highest allowed voltage, $V_{\text{stop}}$ . Inserting $V_{\text{stop}}$ in Eq. (121) yields:

[TABLE]

Thus a single level QD in the sequential tunneling regime always achieves $\eta_{\text{max}}=\eta_{C}$ , regardless of the the temperature difference, $\theta_{0}$ and the distance between $\mu$ and $E_{1}$ bib:linke ; bib:humphrey . Using Eq. (120) we can compute the power when $\eta$ is maximum, i.e. when $V=V_{\text{stop}}$ : this yields $P=0$ . These results agrees with the expectation that Carnot’s efficiency is reached when the heat exchange is “reversible”, thus when the power is vanishingly small.

Let us now study the efficiency at peak power, $\eta(P_{\text{peak}})$ . $P_{\text{peak}}$ is obtained by maximizing the power with respect to $\bar{\Delta}_{\text{min}}$ and $V$ , at fixed $\Delta T>0$ , imposing $\bar{\Delta}_{\text{min}}>0$ and Eq. (123). By imposing this request we obtain two coupled equations that cannot be solved analytically esposito2009 . However, the Fermi function is always evaluated when the argument is positive, and we approximate it with its exponential tail. By doing so, we obtain

[TABLE]

where $\bar{e}$ is Napier’s constant. These equations provide an approximate expression of $P_{\text{peak}}$ and $\eta(P_{\text{peak}})$ for a single level QD, valid for any reservoir temperatures. Note how $\eta(P_{\text{peak}})$ only depends on $\eta_{C}$ , while $P_{\text{peak}}$ depends on both $T_{1}$ and $T_{2}$ through $\eta_{C}$ and $\Delta T$ . Note $\eta(P_{\text{peak}})\rightarrow\eta_{C}/2$ as $\eta_{C}\rightarrow 0$ , as expected from the fact that $ZT\to\infty$ for a narrow transmission probability bib:mahan .

As we can see in Fig. 14, there is a good agreement between $\eta(P_{\text{peak}})$ given in Eq. (125) and a numerical calculation: the analytic expression slightly underestimates $\eta(P_{\text{peak}})$ . Furthermore, we can see that the efficiency at peak power goes beyond $\eta_{CA}$ , while it remains under $\eta_{SS}$ . We have also verified that $P_{\text{peak}}$ given in Eq. (125) is in good agreement with the numerical calculation.

To further assess the validity of the approximate analytical formulas (125), we expand them around $\eta_{C}=0$ and $\eta_{C}=1$ , and compare the obtained results with the exact expansions known in these limiting cases esposito2009 . An expansion of Eq. (125) around $\eta_{C}=0$ yields

[TABLE]

where $T$ is the average temperature in the linear response regime. Our result for $\eta(P_{\text{peak}})$ has to be compared with exact expansion of Ref. esposito2009, ,

[TABLE]

As we can see the first two orders are exactly the same, while the third order only has $\approx 6\%$ correction, confirming that Eq. (125) slightly underestimates the exact result. The expansion of $P_{\text{peak}}$ instead confirms that in the linear response regime the peak power depends on $\Delta T^{2}$ . Furthermore, evaluating the coefficient numerically yields $P_{\text{peak}}\approx 0.105\gamma k_{B}T(\Delta T/T)^{2}$ , which is in very good agreement with the result obtained in Eq. (52) for a multilevel QD in the quantum limit linear response regime (note that $P_{\text{peak}}=Q^{*}\Delta T^{2}/4$ ).

An expansion of Eq. (125) around $\eta_{C}=1$ yields

[TABLE]

to be compared with the exact expansion esposito2009

[TABLE]

There is a good agreement between the two formulas, and we can see that Eq. (125) slightly underestimates $\eta(P_{\text{peak}})$ . We stress that the expression of $P_{\text{peak}}$ for $\eta_{C}\approx 1$ , Eq. (129), shows that the peak power is proportional to $\Delta T$ , as opposed to $\Delta T^{2}$ as in the linear response regime. Furthermore the peak power approaches its maximum value given by $P_{\text{peak}}=\gamma k_{B}\Delta T/(1+\bar{e})$ when $\eta_{C}=1$ .

At last, we have computed numerically the maximum power $P_{\text{max}}$ in the case of a two-level non-interacting QD. $P_{\text{max}}$ is computed by maximizing the power only over $V$ ; the peak of $P_{\text{max}}$ as a function of $\mu$ will thus yield $P_{\text{peak}}$ . In Fig. 15 $P_{\text{max}}/P_{\text{peak}}$ is plotted as a function of $\mu$ for $\Delta E=20\,k_{B}T$ and $\Delta E=2\,k_{B}T$ respectively in panel (a) and (b). Each curve corresponds to a different value of $\eta_{C}$ , starting from the linear-response case (black solid curve) to the extremely non-linear behavior (black dashed curve) at $\eta_{C}=0.99$ . In panel (a), representing the quantum limit, all curves nearly coincide with the linear-response one (apart for very large values of $\eta_{C}$ ). For each value of $\eta_{C}$ the height of the four peaks is equal, contrary to the case where interaction is present (see Fig 11). This is to be expected, since the two energy levels do not influence each other in the quantum limit in the absence of interaction, so the behavior is essentially dictated by a single level. In the bottom panel instead, representing a case away from the quantum limit, the external peaks are much higher with respect to the internal peaks, especially for small values of $\eta_{C}$ . This is due to the fact that $\Delta E$ is of the order of $k_{B}\bar{T}$ and therefore we cannot consider a single energy level at a time, and this produces the peak asymmetry.

Bibliography107

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) M. S. Dresselhaus, G. Chen, M. Y. Tang, R. G. Yang, H. Lee, D. Z. Wang, Z. F. Ren, J. P. Fleurial, and P. Gogna, Adv. Mater. 19 , 1043 (2007) . · doi ↗
2(2) G. J. Snyder and E. S. Toberer, Nat. Mater. 7 , 105 (2008) .
3(3) C. J. Vineis, A. Shakouri, A. Majumdar, and M. G. Kanatzidis, Adv. Mater. 22 , 3970 (2010) .
4(4) G. Benenti, G. Casati, K. Saito, and R. Whitney, ar Xiv:1608.05595 .
5(5) A. I. Hochbaum, R. Chen, R. D. Delgado, W. Liang, E. C. Garnett, M. Najarian, A. Majumdar, and P. Yang, Nature 451 , 163 (2008) .
6(6) L. D. Hicks, T. C. Harman, and M. S. Dresselhaus, Appl. Phys. Lett. 63 , 3230 (1993) .
7(7) G. D. Mahan and J. O. Sofo, Proc. Natl. Acad. Sci. USA 93 , 7436 (1996) .
8(8) C. W. J Beenakker, Phys. Rev. B 44 , 1646 (1991) .

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Thermoelectric properties of an interacting quantum dot-based heat engine

Abstract

pacs:

I Introduction

II Multi-level interacting QD

II.1 Kinetic equations

II.2 Charge, energy, and heat currents

II.3 Detailed balance equations

II.4 Level balance equations

II.5 Output power and efficiency

III Linear response regime

III.1 Quantum limit

III.2 Two-terminal system

III.2.1 Power factor

III.2.2 Figure of merit

III.3 Comparison with a non-interacting QD

III.4 Three-terminal system

IV Non-linear response

IV.1 Non-linear Seebeck and Peltier coefficients

IV.2 Efficiency and output power

V Conclusions

Acknowledgements

Appendix A The kinetic equations always allows a non-trivial solution.

Appendix B The DBEs are not consistent in general

Appendix C Conditions of validity of the linearized DBEs

C.1 Two energy level system (L=2L=2L=2)

C.2 LLL energy level system

C.2.1 Recursive decomposition

C.2.2 The independent part of B(L)B(L)B(L)

C.2.3 Consistency of the linearized DBEs

Appendix D Non-interacting QD: Non-linear effects

C.1 Two energy level system ( $L=2$ )

C.2 $L$ energy level system

C.2.2 The independent part of $B(L)$