Thermoelectrics of Interacting Nanosystems -- Exploiting Superselection instead of Time-Reversal Symmetry

Thermoelectric transport is traditionally analyzed using relations imposed by time-reversal symmetry, ranging from Onsager's results to fluctuation relations in counting statistics. In this paper, we show that a recently discovered duality relation for fermionic systems -- deriving from the fundamental fermion-parity superselection principle of quantum many-particle systems -- provides new insights into thermoelectric transport. Using a master equation, we analyze the stationary charge and heat currents through a weakly coupled, but strongly interacting single-level quantum dot subject to electrical and thermal bias. In linear transport, the fermion-parity duality shows that features of thermoelectric response coefficients are actually dominated by the average and fluctuations of the charge in a dual quantum dot system, governed by attractive instead of repulsive electron-electron interaction. In the nonlinear regime, the duality furthermore relates most transport coefficients to much better understood equilibrium quantities. Finally, we naturally identify the fermion-parity as the part of the Coulomb interaction relevant for both the linear and nonlinear Fourier heat. Altogether, our findings hence reveal that next to time-reversal, the duality imposes equally important symmetry restrictions on thermoelectric transport. As such, it is also expected to simplify computations and clarify the physical understanding for more complex systems than the simplest relevant interacting nanostructure model studied here.


Motivation and general outline
A thorough understanding of the thermoelectric operation of basic circuit elements such as quantum dots is important for future energy-and information-technologies, see, for example, the review articles [1,2] and references therein. Their nonlinear operation, due to large temperature and voltage gradients, is nontrivial and has only recently received more attention [1][2][3][4][5][6][7][8][9][10][11][12][13], the results indicating new possibilities for thermoelectric applications. To get a better grip on this nonequilibrium regime, the implications of time-reversal symmetry -often exploited in thermoelectrics -have been expressed in fluctuation-relations using the powerful tools of counting statistics [14][15][16][17][18]. This is helpful since in strongly confined devices the coupled nonlinear transport of charge and heat is particularly rich and correspondingly difficult to compute. The reason for this is the energy-dependent transmission caused by the system itself, having discrete quantum states and strong Coulomb interaction, rather than by its coupling. In particular, the effect of Coulomb interaction is multifaceted, since local interactions can both store energy and modify the transfer of energy by individual particles. For applications this can be both damaging and helpful: interaction can limit the efficiency since it provides a way to transfer heat without inducing a charge current, thereby leading to leakage; yet, it can also modify the effective level structure in an advantageous way [6,19] or be used for energy harvesting in multi-terminal devices [20][21][22] as shown experimentally [23][24][25].
One reason why the discussion of efficiencies and performance of thermoelectric devices has been mostly confined to the linear regime of operation is the powerful Onsager reciprocity [26][27][28]. The latter is implied by time-reversal symmetry and it allows to make straightforward statements about efficiencies. However, recently a new general symmetry relation has been discovered for electronic systems [29] that is independent of time-reversal symmetry. Its formulation starts from the very general observation that the dynamics of any system obeying linear superposition of its states is characterized by its time-evolution modes, the right eigenvectors of the evolution generator. The amplitude of excitation of these modes is governed by the corresponding left eigenvectors. For closed systems, these mode and amplitude vectors of the Schrödinger dynamics ∂ t |φ = −iH|ψ are the energy eigenkets and bras, respectively, which are trivially related to each other by taking the Hermitian-adjoint (since H = H † ). However, for open systems no such simple general relation exists due to the more complicated nature of nonunitary time evolution. Nevertheless, some of us have shown [29] that for any fermionic open system in the wide-band limit, a quite general duality exists between the modes and amplitudes of the time-evolution kernel. This holds even when the reduced density operator of the open system obeys a kinetic equation with the time-nonlocal form ∂ t ρ(t) = −i[H, ρ(t)] + t 0 dt W (t,t )ρ(t ). Importantly, the duality relies on a fundamental principle other than time-reversal, namely the superselection rule based on fermion parity [30][31][32][33]. This principle is obeyed by all fermionic systems, and therefore the duality is expected to have far-reaching physical consequences.
For thermoelectric transport, the implications of this additional general principle are to date largely unexplored. Fully exploiting the new duality relation, the present paper revisits the thermoelectric response of a weakly coupled, but strongly interacting single-level quantum dot subject to both electrical and thermal biases, both in the linear and nonlinear regime. We derive new results, significantly clarify known results, and thereby -most notably -outline a completely new approach for the analytical computation of thermoelectrics of interacting nanoscale systems. In Refs. [29,34], it was already realized that the duality relation between decay modes and their excitation amplitudes would naturally express itself in the transient charge and heat currents as the system decays to a new stationary state after a sudden switch (quench). This paper shows that the duality is a similarly powerful tool for the experimentally more accessible stationary thermoelectric transport. Indeed, we demonstrate that both in the linear and nonlinear stationary regime -after some judicious choices -an expansion in open-system evolution eigenmodes is advantageous. The compact analytical formulas we obtain in this way offer several new interesting insights into important measurable thermoelectric transport quantities, and provide a truly nonstandard physical view point. Most importantly and as further outlined below, they expose the remarkable fact that a strongly repulsive system may prominently feature behavior reminiscent of strong attractive interaction in the thermoelectric response coefficients. 1 Since this was not recognized as such due to lack of a proper framework, our paper does not merely present trivial or subjective rewritings, even though we re-derive and re-express some known results. Instead, our approach is uniquely dictated by the physically-motivated evolution-mode decomposition, and the outcomes are virtually impossible to "guess" without the systematics and new intuition provided by the duality.
We finally note that the very general origin and formulation [29,36,37] of the duality allows us to extend our method beyond the weakly coupled quantum dot setup treated here. As our study illustrates, the technical simplifications can be such that one can avoid numerical calculations altogether. We therefore expect that the ideas presented here can also simplify the analysis of more complex systems studied most recently [11,[38][39][40]. This includes situations with broken time-reversal symmetry [34], non-diagonal density operators due to, e.g., noncollinear magnetic fields, and also strongly coupled systems exhibiting time non-local effects. 1 This should be clearly distinguished from the effects studied in [35] for systems which from the start have attractive interaction. Version December 7, 2017 submitted to Entropy

of 38
The paper is structured as follows: after outlining the main ideas and results in general terms [Sec. 1.2] and introducing the microscopic model and the weak-coupling master equation [Sec. 2], we review the minimal details of the duality relation that we will exploit for a gated quantum dot coupled to two electronic reservoirs. The mode-decomposition of the charge and energy current formulas [Eq. (22a)-(22b)] are then analyzed in the linear [Sec. 3] and nonlinear transport regimes [Sec. 4]. The manuscript also contains an appendix which is substantial, not because the new derivations are complicated, but because the steps are nonstandard and deserve to be outlined.

Overview of main ideas and results
By discussing a simple, yet relevant model of a nanoscale system, this paper aims to illustrate how the understanding of the role of Coulomb interaction in thermoelectric transport is advanced by the fermion-parity duality. Before going into further details in Sec. 2.3 it is important to appreciate three of its main ideas in general terms: (1) The duality relation maps the eigenmodes of the system of interest to the amplitudes for a different physical system -the dual system. We will refer to the latter as the inverted system because the duality transformation inverts the interaction (as well as other energies), going from repulsive to attractive and vice versa. For quantum dots, this can be easily visualized by inverting the energy landscape in Fig. 1(a) to that of Fig. 1(c), whose details will be discussed later on. This mapping already explains the seemingly strange occurrence of features of attractive interaction in quantities computed for repulsive systems, as first noted in [29,34]. The straightforward interpretation of such puzzling properties is done resorting to the inverted stationary state, which allows to understand the nontrivial dependence on the original system's parameters from the -often simple and well-understood [41] -physics of the attractive dual model as in Fig. 1(c).
(2) Another reason why the duality clarifies interaction effects is that the "essential" correlating parts of the Coulomb interactions, say, between two orbitals i and j with occupation operatorsn i andn j , respectively, is simply given by parity operators, (1 − 2n i )(1 − 2n j ) = (−1)n i +n j . In fact, correlated electron model Hamiltonians are often directly formulated in terms of the operators on the left hand side. The duality reveals that the total parity operator always corresponds to a special eigenmode of open fermion-system dynamics [36,37], and is hence protected. In simple yet relevant situations, one thereby cleanly separates, throughout the entire calculation, the contributions of the Coulomb interaction into an "essential" correlating part and a nontrivial "average" contribution carried by a charge mode. Since Coulomb interaction is an important source of energy dependence and energy storage in quantum dots, thermoelectrics is thus seen to be intimately tied up with fermion parity and the corresponding duality.
(3) Finally, in the context of thermoelectricity, it is important to emphasize that the duality -in the simple form used here -requires energy-independent coupling between system and reservoir (wide-band limit). This does not mean that it is irrelevant to thermoelectric transport, where properly engineered energy-dependence of the coupling can be of interest for the device operation, see e.g. Refs. [20,21,[23][24][25]. Here, the nanoscale system itself provides the strong energy dependence required for thermoelectric effects, both through strong size-quantization and Coulomb interaction. Models of this kind are relevant in many thermoelectric studies [38,[42][43][44][45] and the duality applies to their description, even when the energy-independent coupling is strong and the temperature is low [29], cf. [36,37]. Also, effective energy-dependent couplings as realized in multi-dot systems [43,44] can be treated in terms of the duality relation presented here. Finally, the duality considerations can be extended [46] systematically to account for the energy-dependence of the coupling.
As a guide to the present paper, we now outline how the above general points turn up in the specific insights of our study, most of which remain hidden when approaching the thermoelectric problem in the standard way: In linear response [Sec. 3], we combine the duality with Onsager's reciprocity derived from time-reversal symmetry. The discussion of the linear response coefficients benefits from the combined insights of both Figure 1. (a) Energy landscape of a single-level quantum dot with repulsive interaction U > 0 in contact with a hot and a cold electronic reservoir with an applied bias voltage. (b) Equilibrium occupation number and parity of the single-level quantum dot as a function of the level position [taken with respect to the right reservoir only, or equivalently obtained for the full system at equilibrium, ∆T = V = 0, as indicated by the "eq" subscript]. The occupation n z,eq shows T -broadened steps at the two Coulomb resonances µ = ε and ε +U, respectively, and the parity p z,eq alternates correspondingly. (c) Energy landscape of the inverted system that is dual to the system of interest in (a). (d) Equilibrium occupation number and parity of the single-level quantum dot in its inverted stationary state. In the dual, attractive system, the charge shows a single double-sized step at µ = ε + 1 2 U, with half the temperature broadening T /2 where two electrons enter or leave the dot by two first-order processes sequential in time, and the parity is essentially always even. Parameters for (b) and (d): U = 10T , V = ∆T = 0, and Γ L = Γ R T .
relations. 2 Even more so, in our particular example, we can show how the duality enforces the Onsager relations in linear response, and restricts how these relations break down beyond the linear regime. To achieve this, we first formulate the linear response in a way that is compatible with the duality. Thereby, we find a simple, explicit expression for the average energy carried by electrons, the tight-coupling part of the heat current. Remarkably, it depends on the mean occupation of the dot in the inverted stationary state: this formula exposes the unexpectedly sharp crossover behavior of the thermopower between well-separated resonances [ Fig. 3] as a two-particle resonance in the dual attractive model. Also, the obvious formal similarity between the linear Ohm and Fourier laws, gets an unexpected twist: Whereas it is well-known that the stationary-state fluctuations of the dot occupation govern the electrical conductance G, the duality reveals that the linear Fourier heat coefficient 3 κ = K − ΠSG, for the heat current in the absence of a charge current, J| I=0 = κ∆T , is dominated instead by occupation fluctuations in the inverted stationary state. As such, the ε-dependence of κ is also governed by attraction, exhibiting the same two-particle resonance as the thermopower. 2 Previously, our duality was similarly combined with the independent well-known Iche-duality [47][48][49] based on an electron-hole transformation, yielding new relations for quantum dots in a magnetic field [34]. 3 The Seebeck coefficient S and the Peltier coefficient Π are introduced in Sec. 3. Version December 7, 2017 submitted to Entropy

of 38
The generally more complicated regime of nonlinear response [Sec. 4] can also be addressed with the duality. For the calculation of the nonlinear Seebeck and Fourier coefficient, we show that an evaluation of equilibrium dot observables -both in the original and dual, inverted system -is often sufficient to obtain the nonequilibrium heat currents. This leads to compact analytical expressions and major simplifications in their explicit calculation, while clarifying the underlying physical picture. For example, the nonlinear Fourier heat is essentially the difference between the parity when the quantum dot is in equilibrium with a single lead on the left or on the right. Although it is well-known that the Fourier heat is carried by electron-electron interaction -even for macroscopic devices [50] -the duality pinpoints precisely which part of the interaction is crucial: it is the parity operator that is entirely responsible for the transferred Fourier heat. Finally, the strong difference between Peltier and Seebeck effects in the nonlinear regime, indicating the breakdown of Onsager reciprocity, can be rationalized completely. The nonlinear Peltier coefficient can be decomposed into an equilibrium Seebeck contribution that stems from the heat transport tightly coupled to the charge transport and a parity-mode contribution. Both of these contributions are sensitive to the attractive-interaction physics of the dual system.

Model, master equation, and duality
2.1. Model, assumptions, and notation We are interested in thermoelectric transport through a single-level quantum dot between two 4 reservoirs labeled α = L,R, as sketched in Fig. 1

(a). The full Hamiltonian
H tot =Ĥ +Ĥ leads +Ĥ tun (2) decomposes into three parts. The dot Hamiltonian describes the correlated dynamics of a single spin-degenerate level ε with Coulomb charging energy U for doubly occupying the level. Heren σ =d † σdσ are number operators for dot electrons of spin σ =↑, ↓ with creation (annihilation) operatorsd † σ (d σ ). The quantum dot is coupled to noninteracting electronic reservoirs α described byĤ with operatorsĉ αkσ andĉ † αkσ for electrons with spin σ and orbital quantum number k. The reservoir energies ε αk and tunneling amplitudes t αk in the tunnel couplinĝ are assumed to be spin-independent. We moreover assume that the relevant coupling strengths Γ α (E) = 2π ∑ kσ δ (ε αk − E)|t αk | 2 determining the typical scale of the tunneling rates are energy independent: Γ α (E) ≈ Γ α . This wide-band limit assumption is crucial for the derivation of the fermion-parity duality in the simple form given below in Eq. (13), see also our earlier remarks in Sec. 1.2 under point (3). Finally, the individual grand-canonical leads α are characterized by the electrochemical potential µ α , and the inverse temperature β α = 1/(k B T α ) entering through the particle / hole Fermi distribution functions f ± α (x) = [e ±β α (x−µ α ) + 1] −1 . We let T = T R and µ = µ R denote the global-equilibrium values of the reservoirs and ∆T and V the applied biases through T L = T + ∆T and µ L = µ − eV . 4 Many statements can be generalized in a straightforward manner to multiple reservoirs.
In the next two sections, we set up the calculation of the stationary charge and heat currents, I α and J α , through the quantum dot which can be written in terms of particle and energy currents, I α N and I α E : We will outline how to compute the currents I α Eq. (4)] from lead α into the dot in a way that exploits the duality. From here on, we set e = k B =h = 1.

Master equation and non-equilibrium currents
Following previous work [29,34,51,52], we outline how to compute the particle and energy currents via the reduced density operatorρ of the quantum dot obtained by tracing out the reservoir degrees of freedom. For weak coupling and high temperatures, Γ α T α , the mixed state obeys a Born-Markov master equation: In Eq. (7), operators are written as supervectors, indicated by round kets:ρ = |ρ). In the following we also need to consider the "superhermitian" conjugate of such a supervector |x) which is conveniently written as bra (x|. This denotes a linear function acting on any |y), i.e., another operatorŷ, as follows: (x|y) := Tr [x †ŷ ]. In this notation, the kernel W is a superoperator whose specific matrix elements are given in App. B for completeness, but which are not required in the following [cf. Eq. (15)]. Such a matrix element ( f |W |i) physically relates to the rate for a tunnel-induced transition, |i) → | f ), between two of the dot states The empty state |0 has energy E = 0, the mixture of the spin states | ↑ and | ↓ has energy E = ε, and the doubly occupied state |2 = | ↑↓ has energy E = 2ε +U. Since the system is fully rotationally invariant 5 and because we focus on the particle-and energy currents, we only need the kernel W in a linear subspace spanned by the three trace-normalized operators given in Eq. (8), see the supplement to [29] for details. Importantly, due to the weak coupling, the kernel is the sum of kernels W α that would be obtained if the system was coupled only 6 to reservoir α: In the same approximation, the currents I α N and I α E can be expressed [34,51] in terms of the reservoir-resolved kernel W α , the number operatorN = ∑ σnσ , as well as the HamiltonianĤ of the dot, and the solution of Eq. (7): In writing the second equality in Eqs. (10) and (11), particle and energy conservation have been used. 7 From this point onward, the standard way to obtain the currents in the stationary limit [namely, where the left hand side of Eq. (7) equals zero] seems straightforward, knowing the Fermi Golden-Rule expressions for the rates ( f |W α |i)/( f | f ): (i) Find the non-equilibrium state |z) by solving the stationary limit of Eq. (7), W |z) = ∑ α W α |z) = 0 ("zero mode"), and normalize it to (1|z) = Tr [ẑ] = 1. (ii) Plug |ρ) = |z) into the current 5 The rotational symmetry of the full Hamiltonian (2) causes the master equation for the probabilities in the states (8) to completely decouple from the local spin-dynamics, including the coherences. 6 Regarding the use of the reservoir index α, we note that the particle and energy current in reservoir α depend on the properties of all reservoirs. This is in contrast to, e.g., the kernel W α , whose contributions are due to a single reservoir α only. We emphasize this difference by choosing superscripts for the former case and subscripts for the latter. 7 The second step in the energy current Eq. (11) is only valid in the weak-coupling limit, see remarks in Refs. [29,51] and Ref.
[53]. equations (10) and (11). While this procedure yields explicit expressions for the non-equilibrium currents, it often provides only limited analytical and physical insights, in particular for thermoelectric quantities. Moreover, it unnecessarily complicates the evaluation of bias-derivatives required for the linear-response coefficients, even when using Beenakker's linearization [54,55]. However, inspection of equations (7), (10), and (11) suggests an alternative route based on our recent results [29,34] which we outline next.

Fermion-parity duality and its use in thermoelectrics.
The calculation of the stationary state and the stationary currents are expected to be drastically simplified when expressed in a basis of time-evolution modes, a standard practice when solving linear systems of equations. Looking at Eq. (12), three different eigenvector bases suggest themselves: those of the two separate kernels W α (α = L,R) and of the total kernel W = W L + W R . A key technical point of the present paper is that even though the evaluation of the currents (12) depends on the total kernel W through its zero mode |z), it often turns out that a mode decomposition of the separate kernels W α suffices to compute |z). This provides an important advantage: the kernels W α of the separate leads are by definition always equilibrium kernels. We can thus express and understand nonequilibrium currents in terms of equilibrium quantities which have much simpler dependencies on parameters.
There are three reasons why this decomposition is possible. First, the weak-coupling approximation allows the decomposition given in Eq. (9). Second, our study focuses in most parts either on the linear response regime [Eq. (24)], or on the nonlinear regime but with a vanishing average charge current. We will return to this point at the end of this section. The third, most crucial ingredient is the fermion-parity duality, which strongly restricts the form of the kernels W and W α and thereby identifies the optimal variables in which to express the currents.
We first discuss the duality for a quantum dot coupled to only one lead α, allowing us to take over the expression for W α given in Ref. [29] for this case. There, it was shown that there is a duality relation between the kernel W α and its Hermitian conjugate at inverted energies: Here, W α (−Ĥ, −µ α ) denotes the kernel of the dual system, in which the signs of all energies in the Hamiltonian H, as well as the electrochemical potential µ α in lead α, have been inverted. This is sketched in Fig. 1(a) and (c). The operator (−1)N is the fermion parity of the dot, giving +1 for even and −1 for an odd occupation number, and its appearance reflects the fundamental origin of the duality, Eq. (13). We refer to Refs. [29,51] for an introduction to fermion-parity superselection, the general duality relation and its special form used here, as given in Eq. (13). The duality relation imposes strong restrictions on the fermionic master equation (7). In fact, these restrictions are so strong that they completely determine the rate matrix W α . For illustration, this construction is carried out explicitly in App. A and directly produces the three eigenvalues −γ mα of W α and the corresponding left and right eigenvectors(m α | respectively |m α ) labeled using m = z, c, p: where the zero modeẑ α does not show up because γ zα = 0. Alternatively, one may directly compute these quantities by diagonalizing the matrix W α employing the duality (13) as done in Ref. [29] and obtain: [parity operator] (15) Inspecting this table, one explicitly sees that the mode vectors and amplitude covectors are cross-related through the parity operator (−1)N and the energy inversion (indicated by subscript "i" onẑ iα and n iα ) that appear in the duality (13). The stationary state -given by the equilibrium grandcanonical ensemble |z α ) =ẑ α = e −β α (Ĥ−µ αN ) / Tr e −β α (Ĥ−µ αN ) [Eq. (A4)] in the here considered weak tunnel-coupling regime -has the eigenvalue 0. Hence, duality relates |z α ) to (p α | with eigenvalue −γ pα = −Γ α , where (p α | is the trace with the inverted stationary state × parity. Similarly, the trace operation (z α | is cross-related to the parity operator |p α ). Finally, since the Fermi function with respect to any chemical potential The crucial point made by the table in (15) is that it reveals the natural variables in which the currents (12) are expressed when inserting the mode expansion (14). These include the expectation values of the charge and parity in the stationary state, entering through the factors (m α |z). However, due to the factors (N|m α ) and (H|m α ), additional expectation values with respect to the reservoir-α resolved equilibrium state [Eq. (A4)] appear: explicitly given as a functions of ε,U, µ α , T α in Eqs. (A10) and (A40), and most importantly, expectation values whereẑ iα is the corresponding inverted stationary state of the dual system [Eq. (A7)], In Fig. 1 (b) and (d), we plot the averages (17) and (18) for n zα and p zα in the right lead, i.e., µ α = µ and T α = T : They are simple, stepped functions of the physical parameters whose shapes are straightforwardly rationalized based on the physics of a repulsive and attractive quantum dot in equilibrium with a single lead α, respectively, as explained in the caption to Fig. 1. Inverting the interaction U -going from Eq. (17) to Eq. (18) -qualitatively changes the parameter dependence of these quantities, as shown in the following. We point out that the standard treatment of the problem completely overlooks that expressions (18) are part of the natural variables in which to express the currents (12), and in which to understand the parameter dependence of the transport. The importance of this insight has been demonstrated for the transient time-dependent response of a quantum dot coupled to a single reservoir in [29]. The final issue that we face in extending this to the stationary nonequilibrium thermoelectric currents (12) is that these currents depend on the stationary state |z) of the system coupled to both leads. This state is the zero mode of the total kernel W , for which the duality holds independently, as shown in Ref. [29]: introducing the lump sum of couplings Γ = ∑ α Γ α . Importantly, while equation (20) itself follows from the lead-resolved duality (13) with W = ∑ α W α , its implications for the eigenvectors of W = − ∑ m=c,p γ m |m )(m| and their eigenvalues are nevertheless non-trivial, since W L and W R do not commute. The nonequilibrium stationary state that we need is, in particular, not generally a sum of lead-resolved stationary states with some lead-resolved 8 Section 3 will, however, show that such a simple decomposition does hold when taking the first temperature-bias or voltage-bias derivative [Eq. (24)] at equilibrium. More remarkably, as exploited in Sec. 4, this decomposition continuous to hold even in the nonlinear regime [Eq. (51)] whenever the system is balanced to maintain zero particle current, I α = 0 [Eq. (41)]. Both cases thus allow us to by-pass the computation of |z) and fully exploit the duality (13) for W α for each lead separately in the way outlined above. This procedure directly ties the non-equilibrium thermoelectric transport to equilibrium physics by introducing the optimal variables (17) and (18) from the start.

Charge and energy currents
Combining all above outlined ideas, we insert the eigenmode decomposition W α = − ∑ m=z,c,p γ mα |m α )(m α | into the current equations (12) and use the table in (15) to obtain the central current formulas of the paper: The particle current (22a) is simply the charge relaxation rate × the excess dot charge relative to its equilibrium value with respect to lead α, where the current is measured. The energy current (22b) shows a more interesting nonstandard decomposition: it has a tight-coupling contribution directly proportional to the charge current, and a contribution which is independent of the charge current (22a), and hence associated to the parity mode. The prefactor of the tight-coupling term involves the energy scale Remarkably, the interaction contribution to this energy is not naturally expressed in the stationary state charge n zα , but instead in the inverted stationary charge n iα with respect to lead α. As expected, E α is approximately ε in the vicinity of transitions between zero-and singly occupied state, and approximately ε +U in the vicinity of transitions between singly and doubly occupied state. However, in the crossover regime between these two resonances it is the physics of attractive interaction in n iα that dominates its behavior, cf. Fig. 1.
The overlap (z iα (−1) N |z) ≈ (z iα |z), contributing to the non-tightly coupled term of Eq. (22b), approximately equals a state overlap since the parity in the inverted stationary state is almost always even due to the attractive interaction of the dual model, see Fig. 1(d). This relation will turn out to be useful for the 8 If one considered more than two leads and obtained 3 linearly-independent vectors already from all |z α ), one could formally expand |z) in these vectors. However, this would not be beneficial, since it would still require knowledge of the explicit form of |z), and the weights λ α would be functions of all system parameters, including the temperature and chemical potential of all other reservoirs. understanding of the characteristic features of the nonlinear Peltier coefficient discussed in Sec. 4.2. For explicit numerical results and plots, it is helpful to write fully expressing it in terms of the measurable quantities introduced in Eqs. (16) to (18). See, e.g., Ref. [52] for the derivation of this relation, and appendices A-C for the explicit calculation of all appearing quantities. The current formulas (22) are valid both in the linear and in the nonlinear response regime and form the starting point of the remainder of the paper. The following sections discuss their implications for the stationary thermoelectric response of the quantum dot from the new perspective offered by the duality (13).

Linear response regime
We start with the investigation of the linear thermoelectric response of the quantum dot to voltage and temperature gradients, |V |, |∆T | T . In this case, we can consider the symmetrized charge and heat currents, I = (I L − I R )/2 and J = (J L − J R )/2 which are standardly written in terms of the Onsager matrix, Here, the diagonal elements L 11 = T G and L 22 = T 2 K present the electric and thermal conductances, G = dI/dV | eq and K = dJ/d∆T | eq respectively, whereas the off-diagonal elements L 12 = T 2 ∂ I/∂ ∆T | eq and L 21 = T ∂ J/∂ ∆T | eq characterize the thermoelectric and the electrothermal responses. Here and below, we denote the evaluation of any quantity q in the equilibrium limit by either of the expressions q| eq = q eq := q| V =∆T =0 . The linearization of our general result (22) for the currents around equilibrium is conveniently found from for derivatives with respect to x = µ α , β α , where we have introduced the affinitiesÂ α := β α Ĥ − µ αN . The derivation of this result [App. C.1] only involves the linearization around equilibrium and nothing else, i.e., no detailed balance or other balancing relations used in previous derivations [54,55] and their extensions [11]. Formulating the linear response on the level of the stationary state has the crucial advantage that it can be combined with the duality. Thereby we can circumvent the cumbersome process of taking derivatives of expectation values and subsequently simplifying them. We instead exploit Eq. (24), the mode expansion given in (15) and two simple orthogonality relations involving the equilibrium state of the original and the inverted model, The first is simply the orthogonality of left eigenvector (p | and right stationary vector |z), here applied at equilibrium, which by duality even holds beyond linear response. The second relation follows from Therefore the second relation is found from Using the above equations, we finally obtain the linear response coefficients of the charge and heat currents (22) [see App. C for the details of this non-standard derivation], with equilibrium decay rates γ c,eq = Γ 2 [ f + (ε) + f − (ε +U )] and γ p,eq = Γ = γ p . These results warrant a detailed discussion given in the following sections. However, what is immediately clear is that even in the linear response regime, the natural quantities in which the response coefficients are expressed are not only the coupling asymmetry Γ L Γ R /Γ 2 , the expected energies (ε, µ, U) and the occupation number n z,eq and its fluctuations in the original dot model, In addition, the occupation number n i,eq [cf. Eq. (22c)] and its fluctuations in the equilibrium state of the dual, attractive systemẑ i,eq appear as well. They enter via the overlap of the equilibrium state and dual equilibrium state, (z i,eq |z eq ) = δ n 2 i,eq · δ n 2 z,eq , as proven in App. C.4. In Fig. 2, we plot the ε-dependence of both fluctuations by evaluating the explicit expression of δ n 2 z,eq as a function of ε,U, µ, T , stated in Eq. (A41). These fluctuations show peaks at ε = µ and ε + U = µ for the original model and at ε +U /2 = µ for the dual attractive model. Using that the fluctuations also follow from an ε-derivative of the equilibrium occupation numbers [Eq. (A42)], δ n 2 z,eq = −T ∂ n z,eq /∂ ε and δ n 2 i,eq = T ∂ n i,eq /∂ ε due to Eq. (29), the observed resonances in δ n 2 z,eq and δ n 2 i,eq are readily understood from n z,eq and n i,eq as shown earlier in Fig. 1(b) and (d).

Electric response
For reference, we note that the linear electric conductance is essentially the charge-relaxation rate γ c,eq × the equilibrium charge fluctuations δ n 2 z,eq plotted in Fig. 2. It shows the well-known Coulomb peaks [54] whenever fluctuations of the charge by one electron are energetically allowed, see Fig. 1(b). Starting from the charge current formula (22a), we note that the conductance (30) derives from the charge-mode of the evolution, as physically expected.

Thermo-electric response and Seebeck thermopower
The charge mode also solely determines the response (27b) of the charge current to a small temperature difference ∆T . To tie this to an energy scale, we consider the thermopower or Seebeck coefficient corresponding to the induced thermo-voltage V | I=0 across the two leads in the absence of a charge current, V | I=0 = S∆T (Seebeck effect). The expression E eq := T S has the unit of an energy, related to the voltage required to counterbalance the thermally induced charge current [56,57]. Our linearized result (27b) gives for this energy Before discussing this result, we immediately note that E eq equals the tight-coupling energy (22c) [for n iα | eq = n i,eq with µ α = µ, T α = T ] which by duality naturally appears in the decomposition of the heat current (22b). As we will detail when introducing the Peltier coefficient in Sec. 3.3, this fact indicates that the Onsager relation L 12 = L 21 is obeyed. However, since we purposely refrain from using time-reversal symmetry, it is first of interest to see how, by using only linear response of the state (24) and the mode-amplitude duality (15), the charge current formula (22a) -apparently of very different form from the heat current -produces nonetheless the same energy scale. This follows by noting that for the linear response coefficient L 12 , only equilibrium correlators are relevant 9 , in particular ĤN eq . The latter 10 introduces the characteristic energy E eq : as a consequence of the orthogonality relations (25), the only contribution to ĤN eq comes from the part of the energyĤ that does not couple to the parity (−1)N [see App. C.3]. The additional appearing correlator N 2 eq is responsible for the µ-shift occurring in E eq .
The energy (31) thus naturally appears by mode decomposition of these correlators, showing that the duality is essential for the understanding of the thermopower. As we noted above [Eq. (22)], the interaction-induced part in E eq , respectively S is governed by the behavior of the average occupation n i,eq in the inverted stationary state. Importantly, this energy is not simply the noninteracting part [n i,eq depends nontrivially on U, cf. Eq. (A10)] or a mean-field-like energy which would involve the average n z,eq /2 rather than (2 − n i,eq )/2. In Fig. 3(a), we plot the Seebeck coefficient, Eq. (31), as a function of the energy level position ε − µ. The sawtooth behavior has been known for a long time [55] and it is traditionally understood as follows: Each linear branch of the curve is explained by considering just one of the two Coulomb resonances and ignoring the other. For each branch, the magnitude of the voltage that can be sustained increases linearly 11 with ε − µ. The shape of the whole curve is then understood as a crossover from the particle-dominated transport of one resonance to the hole-dominated transport of the next one when changing ε − µ: by a continuity argument the curve must cross zero between the two resonances to connect the two branches. However, it remains unclear from this line of arguing whether this is a sharp jump -as here in the linear regime -or a smooth crossover -as in the nonlinear regime discussed later. Moreover, close inspection shows that this sharp jump has an anomalous thermal broadening by half the temperature, T /2.
Both issues are readily understood from our result (31) deriving from duality: the mean energy E eq probed by the thermopower depends on the average charge n i,eq of the quantum dot with inverted energies at equilibrium. As Fig. 1(c) and (d) illustrate, the occupation in the inverted stationary state n i,eq abruptly drops from 2 to 0 once ε − µ exceeds −U /2, as is well-known for impurities with attractive Coulomb interaction [41]. The duality thus 9 The linearization of the particle current (22a) in either x = µ α , T α is proportional to terms of the form (N| d dx |z) eq . By Eq. (24) the ∆T -derivative of |z) translates to a temperature derivative of the Boltzmann factor e −(Ĥ−µN)/T which pulls down bothN andĤ, giving correlators N 2 eq and ĤN eq . 10 Ref. [58] also links the thermoelectric conversion efficiency to mixed charge-heat noise, which in essence relates to the time-nonlocal version of the mixed particle-energy correlator appearing here. 11 In the vicinity of a single resonance, the linear behavior of the Seebeck coefficient can simply be understood from the exact solution of a single noninteracting level with transmission probability T (ω) peaked around ω = ε. Then,  explains the unexpectedly sharp jump in the thermopower in Fig. 3. Moreover, since for U T the occupation for the strongly attractive dual model is well-approximated by we see that the anomalous thermal scale of half the reservoir temperature T /2 appears. This has an intuitive interpretation in the dual model: since the attractive interaction favors a single transition involving an electron pair, each electron feels half of the thermal noise 12 . Note that in our calculation for the attractive dual model, two electrons enter sequentially (in time) by two separate processes with rates described by Fermi functions 13 . Their net effect is that the transition becomes allowed at a single value of ε − µ = −U /2. Since our result (31) expresses the thermopower compactly in the natural and well-understood variable n i,eq , we can easily find a simple yet accurate formula for the level positions at which the thermopower achieves a local minimum (ε − ) and maximum (ε + ). This requires finding ε from the extremal condition 12 Note that here electron pairs are not coherently transported: such coherent pair tunneling appears only in processes of order Γ 2 . In the analysis [59] of the rates for such effects a Bose function appears with 2ε +U in its argument, indicating coherent transport of fermion pairs. Interestingly, we note that this rate as function of ε shows the same halving of the thermal noise as discussed in the present paper. 13 The rate matrix (A22) contains no nonzero rate connecting the empty and doubly occupied state.
The strong attractive interaction in the dual system relative to the anomalous lower temperature scale T /2 makes Eq. (32) an excellent approximation to the left hand side of Eq. (33), already for U 5T . This gives Furthermore, the sensitivity of the thermopower to a gate change in the vicinity of these extrema is relevant for applications 14 . We see that in the crossover regime, the negative slope is dominated by the interactions U T , in contrast to the positive slopes 1/T of the two branches associated with isolated resonances. Fig. 3(b) shows that Eqs. (34) and (35)

Peltier coefficient
The Peltier coefficient Π := ∂ J/∂V | eq ∂ I/∂V | eq determines the heat current generated per transferred particle by a small voltage bias |V | T alone (∆T = 0). It thus also defines an energy scale and our heat-current formula (22b) gives This provides a physically different picture of the energy scale (22c) as compared to the thermopower: it is the characteristic energy carried by electrons across a biased junction. This is reflected by how it is obtained: it is the prefactor of the first, tight-coupling part in the heat current (22b) [cf. Eq. (31)]. Of course, Onsager reciprocity dictates L 12 = L 21 , i. e., that the values of the energy scales obtained from the electro-thermal and thermoelectric response are the same, where L 21 = T ∂ J/∂V | eq and Π = −L 21 /L 11 . However, since we explicitly avoid arguing from time-reversal symmetry, it is of interest to see how the reciprocity relation L 21 = L 12 emerges when only using our duality (13). The crucial point in our approach is that the fermion-parity part of the energy current (22b), carrying part of the interaction energy, is to linear order insensitive to the bias V for any value of the remaining parameters: This follows from the linear-response of the stationary stateẑ [Eq. (24)] together with the orthogonalities (25) [see App. C.4]. The result (37) thus implies that in linear response to the voltage bias, the heat current coefficient L 21 stems from the charge mode only, see also Ref. [64]. Since this is equally true for the thermopower in L 12 , we indeed reobtain Onsager's reciprocity relation. We stress in particular that the tight-coupling of the heat-and charge current is not an approximation in our treatment, but follows from duality. In Sec. 4 we will see that the parity contribution to heat current (22b) does play a role for finite bias, pinpointing how the interesting deviations between Π and S emerge beyond linear response.

Thermal response and Fourier heat conductance
Finally, the thermal coefficient L 22 = T 2 ∂ J/∂ ∆T | eq , given by Eq. (27c), acquires contributions from both terms in Eq. (22b). This suggests that there are two types of energy fluctuations that correspond to linear thermal i,eq of the inverted stationary state at the electron-pair resonance of the attractive dual quantum dot system, see Fig. 2 and Fig. 1(d). The smallness of the Fourier coefficient on the scale determined by U /T and Γ stems from the product of occupation number fluctuations δ n 2 i,eq δ n 2 z,eq . Parameters: transport when a thermal gradient is applied. Indeed, the fluctuation of the grand-canonical energyĤ :=Ĥ − µN determining the equilibrium stateẑ eq ∝ e −Ĥ /T decomposes into two corresponding parts: see App. C.4. Thus, the charge (parity) contribution to the linear heat transport coefficient L 22 [Eq. (27c)] are the charge-(parity-) related energy fluctuations × the charge (parity) rate γ c,eq (γ p,eq ).
In thermoelectric applications, one often wants to know the heat transferred from the hot to the cold reservoir when no electrical power is generated. This purely thermal current J| I=0 ≈ κ∆T for a small thermal bias ∆T at zero charge current is characterized by the Fourier heat conductance, κ = ∂ J| I=0 /∂ ∆T | eq . It is given by T 2 κ = L 22 − L 12 L 21 /L 11 , for which we obtain using Eq. (27). This coefficient stems entirely from the parity-mode contribution to the heat current. As in the Drude theory for metals [50], the pure thermal conductance κ is due to the electron-electron interaction. 15 However, we note that the interaction also nontrivially enters into the other term of the heat current that does not contribute to κ. The heat-current decomposition dictated beforehand by duality thus naturally pinpoints which part of the interaction enters the purely thermal Fourier heat conductance -namely the fermion-parity part. We again stress that both a naive perturbative decomposition of the heat current (noninteracting part + interaction corrections) as well as a mean-field decomposition fail to achieve this. The duality furthermore clarifies the parameter-dependence of κ. Being the product of the charge fluctuations in the stationary and in the inverted stationary state, the Fourier coefficient might be expected to reflect the fluctuations δ n 2 z,eq of the actual system, for example in its ε-dependence. However, these fluctuations that we plotted in Fig. 2 are unable to account for the single resonance that κ shows at ε − µ = −U /2 in Fig. 4. Instead, this peak entirely comes from the fluctuations δ n 2 i,eq in the dual model, where due to the attractive interaction only a transition from zero and double dot occupation occurs precisely for this level position, cf. Fig. 5. Interestingly, the anomalous thermal broadening T /2 in the attractive model [Eq. (32) ff.] is crucial to make these fluctuations Figure 5. Sketch of the processes leading to Fourier heat transfer in the thermally broadened region around the particle-hole symmetry point, ε − µ = −U /2, for (a) the normal and (b) the inverted quantum dot models. In both cases shown, the energy U flows from the hot lead to the cold one without net charge transfer from the left to the right lead. dominate: it leads to a difference in the exponential thermal suppressions of the two fluctuations when ε is varied, as shown in the inset of Fig. 2. As a result, when multiplying the two curves in Fig. 2, the peak (tails) of δ n 2 i,eq can overcome (suppress) the tails (peaks) of δ n 2 z,eq and thereby determine the shape of the result for κ in Fig. 4. The very good approximation to the peak height for U T , will be derived later on [Eq. (55)] by exploiting the duality directly in the nonlinear regime and linearizing afterwards. Thus, whereas δ n 2 z determines the electric dissipation of the electric conductance, the fluctuation of the inverted dot population, δ n 2 i dominates the thermal dissipation described by Fourier heat. This is the unexpected twist to the formal analogy between Ohm's law and Fourier law that we announced earlier [Eq. (1)]. It appears only when one examines the origin of their coefficients, Eq. (30) and (39), taking into account the fundamental restrictions imposed by duality.
Finally, we show how the duality simplifies the Fourier heat transfer even on a pictorial level. Fig. 5(a) shows how this energy flow can be pictured in the original, repulsive model as a cycle of single-electron tunneling processes, the cycle being restricted to have zero net charge current, I = 0. Two cycles contribute to the Fourier heat, in which the energy ε (ε +U) is transferred from the cold (hot) lead to the hot (cold) one through the lower (higher) resonance 16 of the dot. In either cycle, an amount of energy U is removed from the hot reservoir and released into the cold one. What remains unclear in this picture is how a single sharp peak at ε − µ = −U /2 with a thermal broadening given by T /2 can emerge due to a delicate cancellation between particle and hole processes which are all off resonant. In Fig. 5(b) we show how (in linear response) the same energy transfer can be understood in the inverted dual model in terms of the single resonance available in an attractive quantum dot which can be either empty or doubly occupied. This process clearly shuts off for |ε + U /2 − µ| T /2, implying the fluctuations δ n 2 i,eq vanish and κ with it. 16 corresponding to the addition energy for the transition between zero and single occupation (between single and double occupation)

Nonlinear regime
In the regime where either of the biases V and ∆T is large enough to invoke a nonlinear current response, Onsager's relations are of no help anymore. Although nonlinear fluctuation relations may provide interesting insights [4,65], this would require additional machinery of counting statistics. The fermion-parity duality Eq. (13), instead, can be exploited in the nonlinear regime without further ado: it allows us to simplify the derivation and improve our understanding of the nonlinear Seebeck, Peltier, and Fourier coefficients obtained from the full non-equilibrium currents (22). We note that since now heat currents are not anymore conserved due to the finite Joule heating, it is only meaningful to consider lead-resolved heat currents J α for α = L or R.

Thermo-electric response and Seebeck thermopower
We start by analyzing the nonlinear thermopower / Seebeck coefficient S nl = V | I=0 /∆T . The required thermo-voltage V | I=0 is obtained by solving I = 0 for V using the currents (6). The form of the charge current (22a) shows that this is equivalent to maintaining equal lead-resolved and, therefore, equilibrium occupations on the dot: This simplification is ultimately a consequence of the fact that for the spin-degenerate quantum dot model considered here, duality (13) dictates the charge mode to be an exact eigenmode both in and out of equilibrium, see (15). Since we consider µ R = µ to be fixed, we can readily solve 17 the balance equation (41) for µ L ≡ µ −V and obtain [App. D] the nonlinear thermopower, using the explicit expressions for n zL and n zR given in Eq. (A10). While the result could be written in terms of Fermi functions directly, a particularly insightful form is where n iR is the ε-dependent inverted stationary charge number with respect to the reference chemical potential µ and temperature T in the right lead. The representation in terms of the dual occupation number is motivated by its connection to the linear limit, in which n i,eq = n iR governs the interaction-related contribution to the Seebeck coefficient for a fundamental reason. This allows us, in the following, to exploit the simple single-step behavior and anomalous thermal scale T /2 of the dual charge [Eq. (32)] in order to further analyze S nl .
Let us first consider small to moderate temperature gradients 0 < ∆T T . Plotting S nl as a function of ε − µ in Fig. 6, we see in the here studied regime U T that the nonlinear thermopower maintains the characteristic ε-dependent shape from the linear regime, exhibiting in particular the qualitative signature of n iR . The only important impact of the increasing thermal bias ∆T is that the slope of S nl at ε + U /2 = µ becomes increasingly less negative, and correspondingly, local extrema are less pronounced in Fig. 6(b). To analytically understand this behavior, we realize that due to ∆T T and due to the sharp step of n iR (ε) around ε − µ = −U /2, we can efficiently carry out a linear expansion of Eq. (42) in ε − µ by first expanding S nl in 1 − n iR around 1 − n iR = −1, 0, 1 for ε − µ close to 0, −U /2, −U respectively: For ε − µ close to 0, −U, the flat ε/U-dependence and the anomalously small thermal broadening of n iR by T /2 [Eq. (32)] lead to a complete suppression of the ∆T -dependent non-equilibrium terms. By contrast, the sizable linear ε/U-dependence of 1 − n iR with slope −U /T around ε − µ = −U /2 [Eq. (35)] does introduce a relevant ∆T -dependence of the slope of S nl linearized in ε, given by the inverse effective temperature As expected, taking the limit ∆T → 0 in Eq. (43) immediately gives back the linear Seebeck coefficient (31). We also note that a useful demarcation of the regime where the second line in Eq. (43) is a good approximation is obtained by finding the crossing points of the linear ε-expansions in all three regimes: The approximations in Eqs. (43) and (45) will be helpful for the understanding of the nonlinear Fourier coefficient in Sec. 4.3.
Let us now address larger temperature gradients ∆T T , approaching and exceeding the interaction strength U. An interesting feature of this regime is that the roots of S nl (ε) close to the Coulomb resonances are shifted towards the root at ε − µ = −U /2, and that all roots can even merge entirely at this single level-position for large enough ∆T . This both experimentally [9] and theoretically [66] studied effect is captured by the analytic expression (42) of the thermovoltage, as clearly visible in Fig. 6(c). Our duality makes explicit that there is no reason that the thermopower roots at ε − µ = 0, −U should remain fixed; only the position of resonance at ε − µ = −U /2 is dictated by the dual model. Indeed, the linear thermopower (31) is not exactly 0 at ε − µ = 0, −U, although this effect is exponentially suppressed with large interaction [0 < (2 − n i,eq )/2 ∼ e −U /T at ε = µ].
For large temperature gradients ∆T /T 1, the effect can, however, clearly be seen: figure 6(d) displays a sizable root shift on the scale U /2 as a function of ∆T /U, obtained numerically from Eq. (42) for different T /U. We also observe that all roots merge at ε − µ = −U /2 at a relatively large temperature gradient, which we define as ∆T 0 /U. To get an approximate analytical description of this quantity, we start from the condition dS nl /dε = 0 at ε − µ = −U /2, which sets the point beyond which only one root can exist. Since this typically happens at large ∆T /T , we can Taylor expand the derivative in T /∆T at fixed T /U to analytically solve for Equation (46) is in good agreement with the numerics shown in Fig. 6(d). Altogether, this confirms that the root shift is observable for achievable temperature gradients as long as the the interaction strength U does not exceed the base temperature T by more than an order of magnitude.

Electro-thermal response and Peltier coefficient
The nonlinear Peltier coefficient Π α nl = (J α /I α )| ∆T =0 determines the heat transferred per charge in response to a finite voltage bias V (∆T = 0). It therefore quantifies one of the desired thermoelectric effects: the ability of the quantum dot system to electrically cool one of the contacts (via a heat current out of this contact). Due to the absence of heat current conservation under nonequilibrium conditions, the nonlinear Peltier coefficient depends on the lead α that it refers to, Π L nl = Π R nl , and furthermore differs substantially from the Seebeck coefficient, as we now discuss.
Since the Peltier coefficient is not restricted by charge-current balancing 18 , we can not express it entirely in lead-resolved equilibrium variables. Nevertheless, we can exploit our duality to a large extent. First, we decompose the Peltier coefficient into a part that stems from the heat current contribution that is tightly coupled to the charge current, Π α nl,tc [first term in Eq. (22b)], and the non-tightly coupled contribution, Π α nl,ntc , stemming from the fermion-parity mode [second term in Eq. (22b)], The first term of Eq. (22b) couples the charge current to E α , the average energy (22c). This is the same energy that determines the linear thermopower (Seebeck coefficient), but with respect to µ α . Therefore, and Π L nl,tc is just the linear thermopower curve shifted by the applied bias V , in contrast to Π R nl,tc , which is not shifted. This can be seen in Fig. 7(a) and (b), where we plot the Peltier coefficient and its decomposition, the dashed curves showing the tight-coupling part. It means that, interestingly, the nonlinear Peltier coefficient, although not related to the nonlinear Seebeck coefficient, is actually related to the linear Seebeck coefficient.
The remaining non-tight-coupling contribution Π α nl,ntc is completely due to the parity mode, Using Eq. (22d), this can be expressed in the average occupation and parity in the stationary state |z) and α-resolved inverted stationary state |z iα ). In order to make further quantitative progress, the evaluation of the expectation values n z and p z with respect to the stationary state |z) cannot be avoided (see App. B). However, we stress that Eq. (49) together with Eq. (22d) already present a significant simplification relative to the standard way of writing the same result. Moreover, we can still further exploit the duality in a qualitative discussion. For example, the results in Fig. 7 show that Π L nl,ntc contributes only in a limited interval of ε values. This can be explained using only the observation that this contribution depends on the overlap of the stationary stateẑ and the α-resolved inverted equilibrium stateẑ iα , as given in Eq. (49). Intuitively, this overlap is a measure of how much the two mixed statesẑ andẑ iα resemble each other, modulo parity signs. Taking µ L = µ −V , the inverted equilibrium statê z iL switches from an empty state for ε − µ ≤ −U /2 −V to a doubly occupied state for ε − µ ≥ −U /2 −V . In contrast, the nonequilibrium stationary stateẑ has a less simple bias dependence and determines the evolution of the panels of Fig. 7(b), where V was chosen to be negative, |V | = −V : (Left): For |V | U /2 the overlap (49) is nonzero only at ε +U /2 = µ up to thermal broadening. (Center): For U /2 < |V | < U, there exists a region ε − µ ∈ [0, |V | −U /2], where the stateẑ is a mixture of an empty and singly-occupied state while the inverted stateẑ iL is empty, leading to a nonzero overlap (49). (Right) For |V | > U there is a double step. Here, the inverted stateẑ iL is empty for all ε − µ ≤ |V | −U /2. The double step in the overlap (49) arises due to the stateẑ being empty with probability p 0 ≈ 1/4 and 1/3 for ε − µ ∈ [0, |V | −U ] and [|V | −U, |V | −U /2], respectively.
The duality thus allows to understand the stepwise contributions of Π L nl,ntc and thereby the details of Fig. 7(b): its effect is to shift the tight-coupling part (determined by the Seebeck coefficient, i.e., the "Onsager part") to lower values for moderate bias |V | < U, and to double the saw-tooth behavior for high bias |V | > U.
Overall, a large bias gives more negative values for Π L nl , restricting the regime of effective cooling of the left reservoir to ε − µ > |V |, for which heat is carried out of the left reservoir (equivalently, the right reservoir for this setting can only be cooled if ε − µ < −U).

Thermal response and Fourier coefficient
Finally, we examine the nonlinear Fourier heat coefficient, κ α nl := J α | I=0 /∆T , for lead α = L,R, the ratio of the heat current and the finite temperature bias in the absence of a charge current. By energy conservation 19 , we have κ L nl ≡ −κ R nl . Our heat current formula, Eq. (22b), immediately shows that the nonlinear Fourier heat current is produced solely by the parity mode contribution: In contrast to the expression for the non-tight-coupling contribution to the Peltier coefficient, Eq. (49), the above expression (50) can be fully analyzed in terms of equilibrium quantities, dictated by the duality. We therefore make use of the remarkable fact that under current balanced conditions (41), the full non-equilibrium stationary state |z) simplifies to a sum of the lead-resolved equilibrium states |z α ): This is shown App. D, in fact for any number of leads α. With this, we find Remarkably, it shows that the nonlinear Fourier heat coefficient can be rationalized in terms of two relatively simple equilibrium observables of the quantum dot, the lead-resolved parities p zL and p zR . As previously pointed out, the parity p zR with respect to the right lead, plotted Fig. 8(c), simply equals the equilibrium parity p z,eq , since µ R = µ and T R = T are kept fixed in the right lead. This equilibrium parity simply changes sign at both resonances, ε − µ = 0, −U. In contrast, the required value p zL | I=0 is obtained by evaluating the explicit expression (A40) of p z,eq = p zR at temperature T + ∆T , and at a chemical potential shifted compared to µ R = µ by the nonlinear thermovoltage V | I=0 = S nl ∆T , which depends nontrivially on the parameters [Eq. (42)]: In the following, more detailed analysis of Eq. (52), we can exploit (53) in the regime of moderate temperature gradients 0 < ∆T < T . Namely, an inspection of Figs. 8(c) and (d) shows that the parity with respect to the left lead, Eq. (53), is expected to take a constant value, equal to the value of p zR = p z,eq at the electron-hole symmetric point but at a different temperature T + ∆T , in a ∆T -dependent range around ε − µ = −U /2. And indeed, explicitly evaluating the thermovoltage (42) in the piecewise-linear approximation (43) and for U ∆T /(T (T + ∆T )) 1, we find With the help of equations (52) to (54), we can now fully understand the plots of κ L nl as a function of the dot level ε as shown in Fig. 8(a) and (b) for different values of the Coulomb interaction. For increasing thermal bias ∆T , the linear response peak at ε − µ = −U /2 [ Fig. 4] increases in both height and width, depending on the ratio U /T . Namely, the broadened peak around ε − µ = −U /2 for U /T = 10 clearly assumes a plateau shape at large 20 values of U /T , such as U /T = 30 in Fig. 8(b). Now, the approximate expressions for the Fourier coefficient in Eq. (54) tell us that there is a difference between the parities p zL and p zR only between the two resonances in the regime delimited by condition (45). The width of this regime, U ∆T T +∆T , first increases linearly with ∆T T and starts to become of order U once ∆T T , as indicated by thin, dashed vertical lines in panels (a) and (b). The difference in parities assumes a constant value in the regime given by U ∆T /(T (T + ∆T )) 1, for which the second line in Eq. (54) is valid. As visible in panel (b) of Fig. 8, this results in the plateau at the maximum value

Conclusion and outlook
We have shown that besides time-reversal symmetry, our duality involving fermion-parity superselection is a crucial general principle for understanding both the linear and nonlinear thermoelectric transport through strongly interacting electronic nanoscale systems. If one is unaware of this, many computed results seem to reflect model-specific features whereas in reality they are fixed from the beginning by the physical restrictions imposed by duality.
Our study illustrates that thermoelectrics is not only interesting for future applications, but also provides new ways in which fundamental effects can be studied: thermal transport properties are particularly susceptible to effects tied to the new duality, since the Coulomb interaction is an additional channel for storing and transporting energy. This naturally brings the fermion-parity operator into the problem, which is an open-system evolution mode that is fundamentally "protected" by the duality.
Concretely, the duality explains why the thermoelectric response of a strongly repulsive system shows features characteristic for attractive interaction. This is particularly visible in the characteristic energy scale of stationary linear transport of heat and charge, which is governed by a resonance occurring at an anomalous energy with unexpected thermal broadening. Both were shown to derive from the occupation number of a quantum dot with inverted Coulomb interaction. In addition, the fluctuations of this dual occupation number were shown to play a similarly important role for the linear heat conductance as the usual occupation fluctuations have for the charge conductance.
The technical advantage offered by the new duality relation is most obvious in the more complicated nonlinear regime. We provided very compact analytical formulas for all transport coefficients, namely the nonlinear thermopower, Peltier coefficient and Fourier heat, allowing for intuitive predictions and analyses of their characteristic features. Strikingly, in most cases, the duality allowed us to express the nonlinear thermoelectric properties in terms of equilibrium variables.
The definite advantages for the analysis and understanding of the thermoelectric response of an interacting single-level quantum dot are expected to extend also to more complex systems with, e.g., multiple levels and contacts. Indeed, the duality applies quite generally beyond the limitations [29] of the present paper, also to low-temperature, strongly-coupled quantum dot systems, which are of current interest [67][68][69] as there are still many interesting thermoelectric properties yet to be explored.
Author Contributions: J. Sp. conceived the idea. J. Sc., A. dM., J. Sp., and J. V. performed the calculations. All authors analyzed the results and wrote the paper.

Conflicts of Interest:
The authors declare no interest in conflicts.

Appendix A. Non-equilibrium master equation kernel from physical principles and symmetries
In this appendix, we derive the non-equilibrium master equation kernel by appealing only to its physicality (dissipativity), the validity of the fermion-parity duality, and the fact that each lead individually is in a local equilibrium state. Moreover, we use that we only consider the first order in the tunnel coupling. This means that only sequential tunneling processes are possible and, due to the full spin symmetry, that the coherent dynamics decouple from the time evolution of the energy eigenstate probabilities. In the following, we can thus restrict the treatment to these probabilities only.

Appendix A.1. Kernel for a single lead
Before we turn to the full, non-equilibrium kernel W , we start by arguing only for any individual kernel W α in the reservoir sum Note that this decomposition can be obtained in the sequential tunneling approximation, in which simultaneous transitions from two reservoirs are neglected. In particular, any W α contains only the dot variables and quantities with respect to the lead α, meaning µ α , β α and Γ α . Now, by construction, each kernel W α must be linear in Γ α . This means that if we turn off the couplings Γ α to all other reservoirs α = α, the complete dynamics of the dot are governed exclusively by W α . Each W α is a physically valid, probability conserving and dissipative kernel in its own right. This implies and furthermore dictates real non-negative transition rates 21 , for transitions from an initial dot energy eigenstate |i) to a different final dot state | f ). Additionally, for reasons that become clear in the next paragraph, we make the -not explicitly required yet plausible 22 -assumption that a sequence of single electron tunneling events to and from lead α eventually connects any possible initial dot energy eigenstate |i) to any other final dot state | f ) = |i): Note now that the first statement in Eq. (A2) means that (1| is a nontrivial left eigenvector of W α to the eigenvalue 0, (1|W α = 0 · (1|. Importantly, the property (A3) then ensures [29] that this eigenvalue 0 is non-degenerate, and that the only trace-normalized right zero-eigenvector |z α ) = 0 is the unique, stationary mixed state of the dot, with truly positive probabilities 0 < P χ < 1 for all dot energy eigenstates |χ). Since the lead itself is assumed to be in equilibrium, and since the grand-canonical ensemble ∼ e −β α (Ĥ−µ αN ) is a possible stationary state of the dot if the latter is in equilibrium with the lead in terms of particle and energy exchange, the unique stationary state |z α ) must be the grand-canonical ensemble: with the partition function Z(ε,U, µ α , β α ) = 1 + 2e −β α (ε−µ α ) + e −β α (2ε+U−2µ α ) . We have used that, due to the spin-degeneracy, the occupation number states and energy eigenstates (8) ofĤ, as well as their dual counterparts, form an orthogonal 23 set of Liouville space vectors spanning the relevant Liouville space: In the wide-band limit, W α belongs to the class for which the fermion-parity duality (13) holds: As in Ref. [29], we assume that neither the physicality (A2) nor the connectiveness of all states (A3) are spoiled by the energy inversion in the dual kernel W α (−Ĥ, −µ α ). In the wide-band limit 24 , this is reasonable because both local energies and electrochemical potential are inverted, and because only energy differences to the electrochemical potential matter for the tunneling rates! The important consequence of this assumption is that the parity rate −Γ α is a non-degenerate eigenvalue of the kernel W α with the left and right eigenvectors (z iα (−1) N | and |(−1) N ), as reported in Eq. (15). Here, we have introduced the inverted or dual stationary 21 Otherwise, it is always possible to find an initial state |ρ 0 ) for which the master equation predicts energy eigenstate projections (i|ρ(t)) < 0 or even Im [(i|ρ(t))] = 0 at some time t, which obviously forbids a probability interpretation. 22 A hypothetical way to violate this condition in the wide-band description of the single-level dot would be to take the zero temperature limit in the Coulomb blockade regime, leading to completely blocked transitions out of the singly occupied state. This is, however, neither physical nor practical, since the zero temperature limit cannot be appropriately described by the Markovian sequential tunneling approximation, and since it is of course impossible to realize T α = 0 exactly in practice. 23 but not orthonormal, since (1|1) = 1/2! 24 For energy-dependent bare couplings Γ(E), this situation would change if Γ(E) had roots on the real axis! state |z iα ) = |z(−H, −µ α )) as defined in Eq. (19) as the stationary state of the dual kernel; it derives from the stationary state |z α ) for lead α by inverting the sign of all local energy scales and of the electrochemical potential. Since the stationary state (A4) is the grand-canonical ensemble, this inverted stationary state can also explicitly be written as By now we have already found two left and right eigenvectors of W α that must be biorthogonal to each other. Since the relevant Liouville subspace and its dual are spanned only by the 3 spin-symmetric physical states written in (8) To efficiently carry out this orthogonalization procedure, we use the fact that also |1), |N − 1), and |(−1) N ) form an orthogonal basis of the Liouville space of interest: This shows that in order to be orthogonal to the parity mode |(−1) N ), (c α | must be a linear combination of (N| and (1|. Enforcing also orthogonality to the stationary state |z α ), applying an analogous procedure for the right vector |c α ), and finally choosing an appropriate normalization factor, such that (c α |c α ) = 1, we obtain Here, n zα = (N|z α ) is the average charge number in the stationary state, and n iα = (N|z iα ) the corresponding average in the dual stationary state: The left vector is thus interpreted as the charge deviation from the stationary average, and the right vector is called charge mode. Since the zero eigenvalue and the parity rate are by construction non-degenerate 25 , both vectors are eigenvectors of W α that can be associated with the decay of the average charge number to its stationary value n zα . The timescale of this decay -the charge rate -is set by the eigenvalue −γ cα to be determined. Collecting all previous results, we arrive at the set of left and right eigenvectors that diagonalizes the kernel, and that is stated in Eq. (15): 25 As noted before, 0 is a non-degenerate eigenvalue of W α . Hence, the vectors (A9) can only be eigenvectors or generalized eigenvectors in a Jordan block of W α to the eigenvalue Γ α , or true eigenvectors to an eigenvalue that differs from 0, −Γ α . Remembering moreover that we have assumed full state connectiveness [Eq. (A3)] and a non-degenerate eigenvalue 0 for the dual kernel W α (−Ĥ, −µ α ), the duality (A6) dictates −Γ α to also be a non-degenerate eigenvalue of W α . Thus, using finally that W α is real, we can conclude that (c α | and |c α ) must be amplitudes and modes of W α to an eigenvalue −γ cα obeying 0 > −γ cα > −Γ α .
with eigenvalues −γ zα = 0, −γ cα , −γ pα = −Γ α . The eigenmode expansion of the kernel reads The explicit expression of γ cα finally follows from the fact that in sequential tunneling, direct transitions between the empty and doubly occupied state are not possible. Using (−1)N |0/2) = |0/2) andN|0) = 0,N|2) = 2|2), this leads to In summary, the equations (A10) to (A13) show that instead of calculating explicit transition rates using, e.g., Fermi's golden rule, an equivalent result for the master equation kernel W α of the single-level dot with one lead can be obtained by appealing to general principles, reasonable assumptions and symmetries. The latter includes in particular the fermion-parity duality.

Appendix A.2. Full multi-lead kernel
The remaining question is what happens in case of a non-equilibrium situation with several leads, at different electrochemical potentials and temperatures. Since the fermion-parity duality holds for the full non-equilibrium kernel W , we can formally obtain its eigenmode expansion by repeating the same procedure as for W α in Sec. A.1, making the same assumptions about the connectiveness and physicality of the (dual) kernel 26 . Similarly to Eq. (A11), this yields the three biorthonormal left and right eigenvectors to the eigenvalues −γ z = 0, −γ c , −γ p = −Γ = − ∑ α Γ α , and where 0 < γ c < Γ restricts the total charge rate γ c . The problem of Eq. (A15) is that it is not yet obvious what exactly the non-equilibrium stationary state |z) and all derived quantities including |z i ) are, since we cannot simply assume the grand-canonical ensemble as for the equilibrium case (A4). Starting from (A1), a straightforward but cumbersome approach to this problem would be to calculate the full kernel explicitly from all single-lead kernels W α , and then explicitly solve W |z) = 0 with (1|z) = 1 to obtain the non-equilibrium stationary state |z). However, our goal here -and in general in this article -is to instead show that for everything we are interested in, we never once need the full expression of |z).
Moreover, |z) must also become the grand-canonical ensemble with respect to the common electrochemical potential µ and inverse temperature β under these conditions, and equivalently for the dual stationary state: 27 with (1|z eq ) = 1 and (1|z i,eq ) = 1. This then finally also implies |c eq ) = |c)| eq (A20) As it turns out, the relations (A20) are almost everything we need to know about the full non-equilibrium stationary state and its dual, apart from the stationarity W |z) = 0 and the orthogonality relation (z i (−1) N |z) = 0 that also hold far away from equilibrium.

Appendix B. Explicit expressions for the kernel
In this appendix, we provide explicit expressions for the full non-equilibrium Kernel W and its stationary state for completeness. We emphasize that these expressions are not needed except for producing plots of the nonlinear Peltier coefficient as shown in Fig. 7, requiring the non-equilibrium parity p.
As mentioned before, the transition rates W f ,i = ( f |W |i)/( f | f ) of the full kernel W can be derived from the sequential tunneling decomposition Eq. (9) of W combined with explicitly known form of the lead resolved kernel W α obtained using duality, see Eq. (A12) combined with Eq. (A7), Eq. (A10) and Eq. (A13). More traditionally, one finds (W f ,i ) α by explicitly using Fermi's Golden rule or first order diagrammatic perturbation theory. In any case, the full matrix of transition rates in the occupation number basis |0), |1), |2) reads 27 While this is intuitively clear, it now follows explicitly from our construction. Namely, we realize from the equations (A10), (A12), and (A13) that at equilibrium, the kernels W α differ only by the coupling prefactor Γ α , such that W α = Γ α W 0 with some common superoperator W 0 . More explicitly, W α | eq = Γ α W 0 (A1) ⇒ W | eq = ∑ α Γ α W 0 = Γ Γα W α | eq . Since W |z) = 0 is a homogeneous equation, and since the stationary state is unique by construction, this implies Eq. (A20).
The stationary state |z) -being the right eigenvector of W to the eigenvalue γ z = 0 is found to be represented by |z) in the basis |0), |1), |2). The inverted stationary state, |z i ), is then simply found by replacing all f + α ↔ f − α . From this expression for the stationary state of the original and the inverted model, the following quantities required for the plot of the nonlinear Peltier coefficient can be obtained from a standard trace operation: n z = Tr{Nẑ}, n i = Tr{Nẑ i }, p z = Tr{(−1)Nẑ}, and p i = Tr{(−1)Nẑ i }.
As intuitively clear, at equilibrium, charge and energy current to every lead α vanish separately: Now the main trick to derive the desired result is to cleverly represent the non-equilibrium stationary state |z) in terms of all local equilibrium states |z α ). Using the complete basis (A11), one can write for any lead α. Multiplying by γ cα , summing over α, and dividing 28 the result by In the second and third step, we have identified the charge current I α N and energy current I α E from lead α as part of the expansion coefficients in order to use the conservation laws (A25) and the equilibrium conditions Eq. (A26) and Eq. (A27) in the following steps.
When taking equilibrium derivatives of |z) with respect to x = µ α or x = β α , we are always free to take the equilibrium limit separately for all terms outside of derivatives (when using the product rule), and to interchange sums with derivatives or the equilibrium limit, since every term in Eq. (A29) is separately continuously differentiable. This finally justifies the following computation: Equation (A30) yields the central relation It reduces equilibrium derivatives of the non-equilibrium stationary state |z) to derivatives of the local equilibrium states |z α ). The latter are in turn straightforward to perform for x = µ α or x = β α . Namely, defining the affinitieŝ A α = β α Ĥ − µ αN , using the exponential form (A4) of |z α ) and the fact that the local dot Hamiltonian conserves the local particle number, [Ĥ,N] = 0, the basic rules of taking derivatives yield The relations (A31) and (A32) together prove the wanted result This is the key ingredient which allows us, in the following steps, to express and relate the linear response coefficients to the additional insights provided by the duality relation. This appeals in particular to the role of the fermion parity and the significance of the dual stationary state |z i ).
Appendix C.2. Linearized charge current and conductance We start with the linear response of the particle current I α N to a variation δ x of the variable x = µ α , T α .
To compute the linear response of the particle current I α N in lead α to a variation of the electrochemical potential δ µ α from the equilibrium, we set x = µ α and realize that only the first term proportional to the average charge fluctuations survives in Eq. (A34). More precisely, we obtain Restricting our considerations to the two-contact case as presented in the main paper the result simplifies to To finally plot the conductance, one needs the explicit expressions for the equilibrium charge fluctuation δ n 2 z,eq as a function of the system parameters, including the equilibrium expectation value ofN 2 . One way to obtain this expectation value is by using the Liouville space identity (A8) and the scalar products 29 This allows us to reduceN 2 tô and to write N 2 eq = (N 2 |z eq ) (A38) = 2(N|z eq ) + 1 2 ((−1) N |z eq ) − (1|z eq ) = 2n z,eq + 1 2 (p z,eq − 1) . (A39) The parity p z,eq is, due to Eq. (A20), calculated by taking the equilibrium limit of the lead resolved parity: With this, we explicitly find δ n 2 z,eq = N 2 eq − (n z,eq ) 2 (A39) = n z,eq (2 − n z,eq ) + 1 2 (p z,eq − 1) where f (x) = f R (x) is the Fermi function with respect to the equilibrium potential µ and temperature T . As stated in the main paper, an alternative way to arrive at the explicit expression Eq. (A41) is to take an ε-derivative of the equilibrium occupation number. Namely, this follows from the relation In summary, the above analysis shows that we can explicitly calculate every ingredient in Eq. (A37): γ c,eq and n z,eq using Eq. (A13) and Eq. (A10) combined with Eq. (A19), and finally δ n 2 z,eq according to Eq. (A41).

Appendix C.3. Charge-Energy correlation and Seebeck effect
We continue by calculating the response of the charge current to a temperature gradient. The first crucial step is to evaluate the equilibrium correlation function NĤ eq by expanding it in terms of the decay eigenmodes and amplitudes (15) taken at equilibrium. Using that all operators, i.e.,N,Ĥ andẑ eq , are hermitian and mutually commute 30 , we can insert the Liouville space identity = (A35) (H|c eq ) · δ n 2 z,eq + (H|(−1) N )(z i,eq (−1) N |N · z eq ).
To further simplify the remaining traces in Eq. (A44), we realize that the scalar products 30 The particle number is locally conserved, [Ĥ,N] = 0 and the equilibrium state is just the Boltzmann factor containingĤ andN in the exponential. and the important orthogonality pointed out in the main text [Eq. (25)] causes the second term of (A44) to vanish. The charge-energy correlator of the dot at equilibrium with the bath is thus simply proportional to the average equilibrium charge fluctuation: The proportionality factor is obtained from the decay mode expansion of the Hamiltonian [29] (H|c eq ) Let us finally calculate the linear response of the particle current to temperature gradients. Evaluating Eq. (A34) with Eq. (A47) and Eq. (A48) gives for With x = T α and dβ α /dT α = −δ αα β 2 α , we furthermore find For the symmetrized charge current response between two contacts to a temperature bias ∆T = T L − T , we then find the relation S = −L 12 /(T L 11 ) = E eq /T used in the main paper: (A51)
Consequently, the equilibrium derivative of the heat current (A52) can be expressed as · (z i,eq |z eq ) for x = µ α , T α . (A56) Let us summarize what we learn from Eq. (A55) and Eq. (A56). The first important message is that for potential gradients, x = µ α , only the first part of the heat current that is tightly coupled to the particle current contributes, since d dµ α (z iα (−1) N |z) eq (A55) = 0.
As claimed in the main text, this ensures compliance with Onsager's reciprocity relation. Namely, the linear response of the heat current with respect to electrochemical potential variations is given by For two contacts, this yields (A59) The final task is to determine the linear response of the heat current to a temperature gradient, involving in particular the second term in Eq. (A56). Motivated by the general fluctuation-dissipation theorem at equilibrium, we first aim to express the state overlap (z i,eq |z eq ) entering the second term in terms of particle number fluctuations. With respect to the stationary and dual stationary state, these fluctuations can be conveniently written as δ n 2 z,eq = (N 2 |z eq ) − [(N|z eq )] 2 (A14) = (c eq |N · z eq ) δ n 2 i,eq = (N 2 |z i,eq ) − [(N|z i,eq )] 2 (A14) = 2(z i,eq · (−1) N · N|c eq ).
Appendix D. Non-equilibrium relations for current-balanced system In this final appendix, we derive the relations used in section 4 to characterize energy transport in a current-balanced non-equilibrium situation, that is, for Since Eq. (A69) means 0 = I α N (22a) = γ cα [n zα − n z ] for all α, and since γ cα > 0 (see appendix A), we have n zα = n z ⇒ n zα = n zα ∀α, α .
Clearly, the converse also holds, i.e., it follows from Eq. (A18) that n z = n zα when assuming n zα = n zα , which leads to I α N = 0 when using Eq. (22a). Let us in the following exploit this equivalence to derive the expansion (51) of the stationary state |z) used in the main text.