---
abstract: |
    We introduce a model for charge and heat transport based on the Landauer-Büttiker scattering approach. The system consists of a chain of $N$ quantum dots, each of them being coupled to a particle reservoir. Additionally, the left and right ends of the chain are coupled to two particle reservoirs. All these reservoirs are independent and can be described by any of the standard physical distributions: Maxwell-Boltzmann, Fermi-Dirac and Bose-Einstein. In the linear response regime, and under some assumptions, we first describe the general transport properties of the system. Then we impose the self-consistency condition, we fix the boundary values $(T_\L,\mu_\L)$ and $(T_\R,\mu_\R)$, and adjust the parameters $(T_i,\mu_i)$, for $i = 1,\dots,N$, so that the net average electric and heat currents into all the intermediate reservoirs vanish. This condition leads to expressions for the temperature and chemical potential profiles along the system, which turn out to be independent of the distribution describing the reservoirs. We also determine the average electric and heat currents flowing through the system and present some numerical results, using random matrix theory, showing that these currents are typically governed by Ohm and Fourier laws.\
    [ 80A20, 81Q50, 81U20, 82C70]{}\
    [ Quantum transport; Quantum dots; Landauer-Büttiker scattering approach; Onsager relations; Entropy production; Random matrix theory; Ohm and Fourier laws]{}
author:
- |
    Philippe A. Jacquet\
    \
    \
bibliography:
- 'refs.bib'
title: |
    [**ThermoElectric Transport Properties\
    of a Chain of Quantum Dots\
    with Self-Consistent Reservoirs**]{}
---

Introduction
============

The study of transport properties of systems out of equilibrium is a fascinating subject in theoretical physics. In particular, various models have been developed to find out what are the underlying microscopic mechanisms giving rise to macroscopic behaviours such as Ohm and Fourier laws [@Jackson; @BLR; @LLP].

One interesting class of models that has been introduced with this purpose is related to the Lorentz gas [@Lorentz; @LS; @LebSpohn; @Wagner; @KlagesPRL; @Rateitschak; @Carlos; @LLM; @EY]. In particular, a classical model for particle and energy transport has recently been investigated by Eckmann and Young [@EY]. In this EY-model, the system is a linear chain of *chaotic* cells, each one containing a fixed freely rotating disc at its centre (see  \[The Classical System\]). Both ends of the chain are coupled to classical particle reservoirs assumed to be in thermal equilibrium at temperatures $T_\L$ and $T_\R$, respectively. Particles are injected into the system from these reservoirs at characteristic rates $\gamma_\L$ and $\gamma_\R$, respectively. Once in the system, the non-interacting particles move freely and collide elastically with the boundaries of the cells as well as with the discs. The particles can also leave the system through these reservoirs.

Note that in this model the discs play the role of energy tanks, allowing the redistribution of energy among the particles in a given cell and thus permitting the system to reach at large times a stationary state satisfying local thermal equilibrium (LTE). As a consequence, this model admits a well-defined notion of local temperature.

One of the main results states that, under some assumptions, the average stationary particle and energy currents through the system are governed by Fick and Fourier laws.

In this paper, we will always assume that the considered system has reached its (unique) *stationary* state and discuss only its *average* transport properties. In particular, we shall only consider average stationary currents and consequently simply call them currents.

In the present work, we shall construct a similar type of model using the Landauer-Büttiker scattering approach [@Landauer; @Landauer1970; @ButtikerFourT; @Buttiker-Scattering; @ButtikerScattering]. This will permit, in particular, to establish an effective *quantum* version of the EY-model.

Basically, in the scattering approach to charge and heat transport, one expresses the (average) electric and heat currents through a conductor in terms of the scattering data of this conductor. Therefore, starting from the EY-model, we shall make the following crucial modifications:

-   The main modification concerns the scattering version of the disc. There are two features of the EY-model that should be pointed out: (i) The exact position of the disc within a given cell is not important, so one may consider the cell represented in  \[The two Systems\] (left). (ii) Once the system has reached a stationary state satisfying LTE, the boxes containing the discs are at equilibrium and consequently the net (average) particle and energy currents into these boxes vanish [@Wagner; @KlagesPRL; @Rateitschak].

    Therefore, in the scattering approach, one can obtain an *effective* description of the discs as follows[^1]: We replace the boxes containing the discs by independent particle reservoirs with characteristic temperatures and chemical potentials such that the net (average) electric and heat currents into these particle reservoirs vanish ( \[The two Systems\] and \[The Quantum System\]). In this manner, the particles in a given cell will be in equilibrium with their respective reservoir, the system will be in local thermal equilibrium.

    Although the boxes containing the discs have finite energy, while the reservoirs have infinite energy, we believe that the transport properties of the system *in the stationary state* (assuming there is one) will be similar in both cases.

    Usually, the condition of zero net (average) particle and heat currents, between the system and a reservoir, is referred to as the self-consistency condition [@BolRV; @RichViss; @Visscher; @Davies; @BLebL; @roy]. If one is only interested in the electric current, the self-consistent reservoirs correspond to the *voltage probes* used in mesoscopic physics, originally introduced to model the inelastic scattering occurring in a conductor [@Buttiker-Reservoir; @ButtikerQuantum; @ButtikerBarriers; @Amato; @ButtikerQuantum2; @Buttiker-Review; @Ando; @Pilgram; @Heidi].

    A *crucial* distinction appears here: In voltage probes *only* the electric currents are set to zero while in this paper the self-consistent reservoirs require to set *both* the electric and heat currents to zero. To our knowledge, this fundamental difference has never been discussed in the formalism of Landauer-Büttiker and consequently makes one of the novel aspects of the present work. In particular, our results generalise works by, for example, Büttiker [@Buttiker-Reservoir; @ButtikerQuantum; @ButtikerBarriers], D’Amato and Pastawski [@Amato], Ando [@Ando] and very recently Roy and Dhar [@roy].

    If the couplings between the system and the self-consistent reservoirs are sufficiently small, then one may interpret the self-consistent reservoirs as ideal potentiometers and thermometers [@Engquist; @Sivan-Imry].

-   The transport properties of the cells will be given in terms of $N$ scattering matrices $S^{(1)},\dots
      ,S^{(N)}$. However, we will work mostly in terms of the scattering matrix $S$ associated to the global multi-terminal system, which is obtained by composing the $N$ local scattering matrices together.

-   In order to have some generality, we will describe the particle reservoirs in terms of the three standard physical distributions: the classical Maxwell-Boltzmann distribution and the quantum Fermi-Dirac and Bose-Einstein distributions. We assume that the particles may carry some charge $e$, so that we may speak of electric currents instead of particle currents. Nevertheless, we shall assume, as in the EY-model, that the particles do not interact with each other.

For obvious reasons, we now call the cells [*quantum dots*]{}. In the EY-model, the outer boundaries of the cells were chosen such that the system was chaotic. We shall at some point impose this property by assuming that the quantum dots are classically chaotic.

The organisation of this paper is as follows. In  \[The model\], we establish the framework of multi-terminal systems in mesoscopic physics, present our model and outline the strategy adopted to obtain the desired chain of quantum dots. Starting with  \[General properties\], we will work in the linear response regime and under the following main assumption: the scattering matrix $S$ associated to the system does not depend on the energy. As we shall see, this assumption, which is a good approximation in certain limiting cases, will lead to some interesting consequences.

More precisely, the main results of this paper are divided into three parts. In  \[General properties\], we present the general transport properties of the system, and discuss in particular the Onsager relations and the entropy production. In  \[The null currents condition\], we solve the self-consistency condition, we fix the boundary values $(T_\L,\mu_\L)$ and $(T_\R,\mu_\R)$, and adjust the parameters $(T_i,\mu_i)$, for $i = 1,\dots,N$, so that the net average electric and heat currents into all the intermediate reservoirs vanish. This condition leads to expressions for the temperature and chemical potential profiles along the system, which turn out to be independent of the nature of the particles, independent of the distribution describing the reservoirs. From these profiles, we work out the average electric and heat currents flowing through the system. Let us emphasise that all these results are obtained rigorously and actually hold for *any* geometry (see  \[The General System\]). These results constitute the main analytical results of this paper.

In  \[Numerical analysis\], we restrict our attention to the linear geometry represented in  \[The Quantum System\]. In analogy with the EY-model, we assume that the quantum dots are classically chaotic and consequently use random matrix theory (RMT) to describe the transport properties of the quantum dots. While we cannot prove the validity of Ohm and Fourier laws in this situation, some approximate derivations are feasible in some particular cases (see  \[Remark Amato\]). We thus turn to numerical simulations and found that Ohm and Fourier laws typically hold (in the sense of RMT) in a chain of quantum dots with self-consistent reservoirs.

The Model {#The model}
=========

We consider a chain of $N$ connected quantum dots, where each dot is coupled to a particle reservoir at temperature $T_i$ and chemical potential $\mu_i$, with $i = 1, \dots, N$. Additionally, the left and right ends of the chain are coupled to particle reservoirs with parameters ($T_\L$,$\mu_\L$) and ($T_\R$,$\mu_\R$), respectively ( \[The Quantum System\]). All these reservoirs, which we also call *terminals*, are independent and inject particles into the system according to some distribution function. We assume that they can also absorb particles without changing their state. In this paper, we shall consider the following three cases ($i \in \{\L, \R, 1, \dots, N\}$): $$\label{Maxwell-Boltzmann}
f^{\rm MB}_i(E) \equiv f^{\rm MB}(E; T_i, \mu_i) = \exp\left(-\frac{E-\mu_i}{k_\B
    T_i}\right)~,$$ $$\label{Fermi-Dirac}
f^{\rm FD}_i(E) \equiv f^{\rm FD}(E; T_i, \mu_i) = \left[\exp\left(\frac{E-\mu_i}{k_\B
    T_i}\right) + 1\right]^{-1}~,$$

$$\label{Bose-Einstein}
f^{\rm BE}_i(E) \equiv f^{\rm BE}(E; T_i, \mu_i) = \left[\exp\left(\frac{E-\mu_i}{k_\B
    T_i}\right) - 1\right]^{-1}~.$$

Here, $E \in [0,\infty)$ is the energy, $k_\B$ is Boltzmann’s constant and $T_i > 0$ is the temperature of the $i$-th terminal. For the chemical potentials there is one subtlety: for the Maxwell-Boltzmann and Fermi-Dirac functions, one has $\mu_i \in (-\infty, \infty)$, while for the Bose-Einstein function, one has $\mu_i \in (-\infty, 0)$. Let us emphasise that, for each $i$, the parameters $T_i$ and $\mu_i$ are *independent* of each other.

We consider that the transport properties of the $k$-th quantum dot are described by a scattering matrix $S^{(k)}(E)$ at energy $E > 0$: $$S^{(k)}(E) = \left(S^{(k)}_{ij;mn}(E)\right)~, \hspace{5mm} k \in \{1, \dots, N\}~,$$ where $i,j$ denote the three possible entrances of the $k$-th dot and the indices $m,n$ denote their corresponding channels. In  \[Numerical analysis\], we shall use these scattering matrices to build the scattering matrix $S$ associated to the *global* system made of $N$ quantum dots. For the present, let us assume that we are given the scattering matrix of the global system: $$S(E) = \left(S_{ij;mn}(E)\right)~,$$ where $i,j \in \{\L, \R, 1, \dots, N\}$ and for each couple $(i,j)$ the indices $m \in \{1, \dots, M_i\}$ and $n \in \{1, \dots, M_j\}$ number the channels in terminal $i$ and $j$, respectively. Therefore, the matrix element $S_{ij;mn}(E)$ is the probability *amplitude* that a particle with energy $E$ incident in channel $n$ in terminal $j$ is transmitted into channel $m$ in terminal $i$.

Being given $S(E)$, one can define the *total* transmission probability $t_{ij}(E)$ that a particle with energy $E$ goes from terminal $j$ to terminal $i$ [@Buttiker-Symmetry]: $$\label{transmission prob}
t_{ij}(E) = \sum_{m=1}^{M_i}
\sum_{n=1}^{M_j} |S_{ij;mn}(E)|^2~.$$

Throughout this paper, we always assume that the following holds.

For each $E > 0$, the complex matrix $S(E)$ is *unitary*.\[Assumption S\]

On the other hand, we do *not* assume the scattering matrix $S(E)$ to be symmetric. Using the unitarity of $S(E)$ one easily obtains the following results: $$\sum_{i} t_{ij}(E) = M_j~, \hspace{2mm} \forall j \hspace{7mm} \mbox{and} \hspace{7mm} \sum_{j} t_{ij}(E) =
M_i~, \hspace{2mm} \forall i~.\label{Property Transmission}$$

[The sums appearing in are over the set $\{\L, \R, 1, \dots, N\}$ of *all* terminals. Unless stated, in what follows every sum will be understood over this set.]{.nodecor}

The Currents
------------

In this paper, we shall always assume that the following holds.

-   The particles do not interact with each other.

-   The system admits a unique stationary state (out of equilibrium).

\[Assumptions\]

Under these assumptions, one can derive for any of the distributions – the expressions for the (average) *electric currents* in a multi-terminal system [@ButtikerScattering; @Buttiker-Review; @Datta; @CA]. One finds $$I_i = \frac{e}{h} \sum_{j} \int_{0}^{\infty} \left[t_{ji}(E) f_i(E) - t_{ij}(E) f_j(E)\right] dE~,\label{Electrical Current 1}$$ where $e > 0$ is the charge carried by the particles and $h$ is Planck’s constant.

[To our knowledge, the expression for the average electric current in a multi-terminal system was first introduced, through formal arguments, by Büttiker in [@ButtikerScattering] and was very recently given a rigorous and general derivation (together with ) in a work by Aschbacher *et al.* [@CA] using a $C^{*}$ algebraic approach.]{.nodecor}

The expression can be interpreted as follows: $t_{ji}(E) f_i(E)$ is the average number of particles with energy $E$ that are transmitted from terminal $i$ to terminal $j$, and $t_{ij}(E) f_j(E)$ is the same but from terminal $j$ to terminal $i$. Therefore, $I_i$ is the *net* average electric current at terminal $i$, counted positively from the terminal to the system.

[To describe the transport of neutral particles (such as neutrons or phonons), it suffices to set $e = 1$ in and one will obtain the particle current.]{.nodecor}

It turns out that most relations that we shall encounter will have the same form in the three considered cases (an example being given by the expression ). Therefore, from now on, all relations not wearing the superscript MB, FD or BE will be assumed to hold in all three cases.

We shall also obtain some expressions which do not depend at all on the distribution function describing the reservoirs. In order to emphasise the fact that such relations do not depend on the nature of the particles, we shall say that they are *universal*. We hope the reader will not be confused with the notion of universality used to refer to properties that are model independent.

Under the Assumption A\[Assumptions\], one can also derive the expression for the (average) *heat current* at terminal $i$ [@Butcher; @CA]. Basically, the idea is to write the particle current, which is given by with $e=1$, as well as the energy current, which is also given by with $e=1$ but with an extra $E$ in the integrand, and then to invoke the first law of thermodynamics $\delta Q_i = dE_i - \mu_i
  dN_i$. One finds $$\label{Current Heat}
J_i = \frac{1}{h} \sum_{j} \int_{0}^{\infty} \left[t_{ji}(E) f_i(E) - t_{ij}(E) f_j(E)\right] (E-\mu_i) dE~.$$

[As we shall see in the linear response analysis, if one considers $\delta \mu_i / e$ and $\delta T_i / T$ as the thermodynamic forces, then $I_i$ and $J_i$ turn out to be the right currents to obtain the Onsager relations.]{.nodecor}

Although the expressions – for the electric and heat currents have a direct physical interpretation in terms of $t_{ij}(E)$, we prefer to work in terms of the quantities $$\Gamma_{ij}(E) = M_i
\delta_{ij} - t_{ij}(E)~,$$ where $\delta_{ij}$ denotes the Kronecker delta. For this, we use the relations to rewrite the integrand in as follows: $$\begin{aligned}
\sum_{j} \left[f_i(E) t_{ji}(E) - f_j(E) t_{ij}(E)\right] &=& f_i(E) M_i - \sum_j f_j(E) t_{ij}(E) \nonumber\\
&=& \sum_j f_j(E) \underbrace{\left(M_i \delta_{ij} - t_{ij}(E)\right)}_{= \
  \Gamma_{ij}(E)}~.\end{aligned}$$ Hence, one can rewrite the currents in the following form ($i\in \{\L, \R, 1,\dots,N\}$): $$\begin{aligned}
I_i &=& \frac{e}{h} \sum_{j} \int_{0}^{\infty} f_j(E)
\Gamma_{ij}(E) \, dE~,\label{Electric Current Def}\\
J_i &=& \frac{1}{h} \sum_{j} \int_{0}^{\infty} f_j(E)
\Gamma_{ij}(E) (E - \mu_i) \, dE~.\label{Heat Current Def}\end{aligned}$$

From one immediately obtains the following properties: $$\sum_{i} \Gamma_{ij}(E) = 0~, \hspace{2mm} \forall j \hspace{7mm} \mbox{and} \hspace{7mm} \sum_{j} \Gamma_{ij}(E) =
0~, \hspace{2mm} \forall i~.\label{Property Conductance}$$

[Note that $\Gamma_{ij}(E) \not = \Gamma_{ji}(E)$ in general.]{.nodecor}

Using the relations one immediately sees that the expressions – for the currents satisfy the conservation of charge and energy, respectively, $$\label{conservation currents}
\sum_{i} I_{i} = 0 \hspace{7mm} \mbox{and} \hspace{7mm} \sum_{i} \left(J_{i}
+ \frac{\mu_i}{e} I_{i}\right) = 0~.$$

[The second relation in emphasises the fact that, in general, energy is conserved but not heat, in the sense that $\sum_{i} J_{i} \not = 0$ in general. However, as we shall see in the next section (see  \[remark heat\]), heat is conserved in the linear response approximation.]{.nodecor}

The Strategy {#The strategy}
------------

In  \[The null currents condition\], we will fix the values $(T_\L,\mu_\L)$ and $(T_\R,\mu_\R)$, and determine the parameters $(T_i,\mu_i)$, for $i = 1\dots,N$, for which the *self-consistency condition* is satisfied: $$I_i = 0 \quad \mbox{and} \quad J_i = 0 \quad \mbox{for} \quad i =
1,\dots,N~.\label{Condition Zero Currents}$$ Actually, this condition will be solved only in the linear response approximation. Then we will work out the remaining electric and heat currents: $I_\L = - I_\R$ and $J_\L = - J_\R$ (these equalities are satisfied in the linear response regime; see  \[remark heat\]). Observe that in this situation, charge and heat dissipation occurs mostly at the left and right boundaries of the system. Indeed, due to the self-consistency condition, the dissipation occurs at first order at the boundaries while at most as second order effects along the system.

Up to that point, the system is described by any scattering matrix $S$ and has therefore no specific geometry (see  \[The General System\]). In order to obtain a system with the *linear* geometry represented in  \[The Quantum System\], and to interpret $I_\L$ and $J_\L$ as the currents flowing through the system, from left to right, we will consider that the *global* scattering matrix $S$ is given in terms of the *local* (linearly ordered) scattering matrices, $S^{(1)},\dots,S^{(N)}$, associated to the $N$ quantum dots.

In  \[Numerical analysis\], we will study numerically the properties of the *composite* scattering matrix $S$, using random matrix theory (RMT), and show that the currents through the chain of quantum dots with self-consistent reservoirs are typically (in the sense of RMT) governed by Ohm and Fourier laws.

General Properties {#General properties}
==================

In this section, we consider that the system is coupled to $N+2$ terminals (we use $\L$ and $\R$ to be consistent with the next sections) and suppose that its transport properties are given in terms of some scattering matrix $S$ (see  \[The General System\]).

In the linear response regime and under the assumption that the scattering matrix $S$ does not depend on the energy, we will determine the general transport properties of the system. All these properties are expected to hold on physical grounds and actually some of them have recently been obtained in full generality [@CA]. Nevertheless, in this paper, we are only interested in some particular situations in which the derivations of these properties are somehow simpler.

Linear Transport Analysis {#Linear transport analysis}
-------------------------

In the remainder of this paper we shall always rely on the following.

We assume that, for $i \in \{\L, \R, 1, \dots,
N\}$, $$\begin{aligned}
T_i &=& T + \delta T_i~,\\
\mu_i &=& \mu + \delta \mu_i~,\end{aligned}$$ where $T > 0$ and $\mu \in \mathbb{R}$ are some reference values, while $\delta T_i \in \mathbb{R}$ and $\delta \mu_i \in \mathbb{R}$ are some “small” perturbations. For definiteness, we assume that all subsequent expressions are first order expansions in $\delta T_i$ and $\delta \mu_i$.\[Assumption Linear transport\]

To study the transport, we consider the first order expansion of any distribution function: $$g(E; T_j, \mu_j) = g(E; T, \mu) + \frac{\partial g}{\partial T_j}(E;
T, \mu) \ \delta T_j + \frac{\partial g}{\partial \mu_j}(E;
T, \mu)  \ \delta \mu_j~.$$ Note that the three considered distribution functions – are of the form: $$\label{form g}
g(E; T_j, \mu_j) = g\left(\frac{E-\mu_j}{k_\B
  T_j}\right)~.$$ Hence, one can rewrite the derivatives of $g$, with respect to $T_j$ and $\mu_j$, in terms of the derivative of $g$ with respect to $E$, and obtain the following result: $$\label{First Order Exp}
f_j(E) \equiv f(E; T_j, \mu_j) = f(E; T, \mu) - \frac{\partial f}{\partial E}(E;
T, \mu) \ \left[\left(\frac{E-\mu}{T}\right) \delta T_j + \delta \mu_j
  \right]~.$$ Substituting this first order expansion into the expressions – for the currents, and using the relations , one sees that all terms containing $f(E; T, \mu)$ vanish and one is left with $$\begin{aligned}
I_i &=& \sum_j L_{ij}^{(0)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(1)}
\, \frac{\delta T_j}{T}~,\label{The electric current}\\
J_i &=& \sum_j L_{ij}^{(1)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(2)}
\, \frac{\delta T_j}{T}~,\label{The heat
  current}\end{aligned}$$ where $$\begin{aligned}
L_{ij}^{(0)} &=& -\frac{e^2}{h} \int_{0}^{\infty} \frac{\partial
  f}{\partial E}(E; T, \mu) \Gamma_{ij}(E) \, dE~,\label{transp 1}\\
L_{ij}^{(1)} &=& -\frac{e}{h} k_\B T \int_{0}^{\infty} \left(\frac{E-\mu}{k_\B T}\right) \frac{\partial
  f}{\partial E}(E; T, \mu) \Gamma_{ij}(E) \, dE~,\label{transp 2}\\
L_{ij}^{(2)} &=& -\frac{(k_\B T)^2}{h} \int_{0}^{\infty} \left(\frac{E-\mu}{k_\B T}\right)^2 \frac{\partial
  f}{\partial E}(E; T, \mu) \Gamma_{ij}(E)  \, dE~.\label{transp 3}\end{aligned}$$

\[remark heat\] [Using the relations , one sees that the first order approximations for the heat currents satisfy $$\sum_{i} J_{i} = 0~.\label{cons heat}$$ In  \[The null currents condition\], we will impose the self-consistency condition, $I_i = J_i = 0$ for $i =
1,\dots,N$, and the conservation laws and will permit to interpret $I_\L = -I_\R$ and $J_\L = -J_\R$ as the electric and heat currents flowing through the system (from left to right).]{.nodecor}

The Transport Matrix {#The transport matrix}
--------------------

Here we establish some basic properties of the transport matrix. From the expressions – for the currents we are naturally led to the following definition.

\[def transport matrix\] We denote by $L$ the transport matrix: $$L = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}L^{(0)} & L^{(1)}\\ L^{(1)} & L^{(2)}\\ \end{array}\hspace{-2mm}\right)}~, \hspace{3mm} \mbox{where}  \hspace{3mm}  L^{(k)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}L^{(k)}_{\L \L} & \L^{(k)}_{\L 1} & \dots & L^{(k)}_{\ LN} &
  \L^{(k)}_{\L \R}\\ L^{(k)}_{1 \L} & \L^{(k)}_{1 1} & \dots &
  L^{(k)}_{1 N} & L^{(k)}_{1 \R}\\ \vdots & \vdots &  \vdots & \vdots
  & \vdots \\ L^{(k)}_{N \L} &
  L^{(k)}_{N 1} & \dots & L^{(k)}_{N N} & L^{(k)}_{N \R} \\ L^{(k)}_{\R \L} &
  L^{(k)}_{\R 1} & \dots & L^{(k)}_{\R N} & L^{(k)}_{\R \R} \\ \end{array}\hspace{-2mm}\right)}~.$$

A quick glance at the expressions – for the transport coefficients $L_{ij}^{(k)}$, remembering that $\Gamma_{ij}(E) =
M_i \delta_{ij} - t_{ij}(E)$, shows that the Onsager relations hold, if $S(E)$ is symmetric (for each $E > 0$), then $L$ is symmetric. More generally, if the system satisfies the micro-reversibility property $S_{ij,mn}(E,B) = S_{ji,nm}(E,-B)$, where $B$ is some applied magnetic field, then $L_{ij}(B) = L_{ji}(-B)$.

In general, the scattering matrix $S$ depends on the energy. Nevertheless, in this paper, we shall restrict our attention to the energy-independent situations. The assumption that $S$ does not depend on the energy is a good approximation in some limiting cases (looking at – one sees that for electrons at low temperature one may consider $S(E) \equiv S(E_\F)$, where $E_\F$ is the Fermi energy, since $(-\partial f/\partial E)(E) \approx
\delta(E-E_\F)$) and leads to interesting consequences. In particular, we will see that the relations among the transport coefficients, $L_{ij}^{(0)}$, $L_{ij}^{(1)}$ and $L_{ij}^{(2)}$, will be very simple and will lead to some universal transport properties. For definiteness, we assume for the rest of this paper that the following holds.

The scattering matrix $S$ does not depend on the energy.\[Assumption Gamma\]

In order to simplify some derivations, without restricting our results too much, it is convenient to assume the following.

We assume that $t_{ij} \not = 0$ for all $i,j \in \{\L,\R, 1,\dots, N\}$.\[Assumption Gamma2\]

[Recalling the properties , one sees that the Assumption A\[Assumption Gamma2\] implies that $t_{ij}
  \not = \min\{M_i,
  M_j\}$, for all $i,j \in \{\L,\R, 1,\dots, N\}$, and this permits to avoid total back-scattering ($t_{ii} \not = M_i$ for all $i$). Although not optimal, the Assumption A\[Assumption Gamma2\] also insures the existence and unicity of the solution of the self-consistently condition .]{.nodecor}

Under the assumption that the scattering matrix $S$ does not depend on the energy, one sees in the expressions for the transport coefficients – that some integrals of the following form appear (with $n = 0, 1, 2$): $$C(n) = - \int_{0}^{\infty} \left(\frac{E-\mu}{k_\B T}\right)^n \frac{\partial
  f}{\partial E}(E; T, \mu) \, dE~.$$ More explicitly, $$C^{\rm MB}(n) = \int_{x_0}^{\infty}
  x^n e^{-x} \ dx \hspace{5mm} \mbox{and} \hspace{5mm}C^{\rm \pm}(n) =  \int_{x_0}^{\infty}
  \frac{x^n e^x}{(e^x \pm 1)^2} \ dx~,$$ where the signs $+$ and $-$ correspond to the Fermi-Dirac and Bose-Einstein cases, respectively, and $x_0 = -\mu/(k_\B T)$ is some reference parameter. Here, $x_0
\in \mathbb{R}$ in MB and FD, but $x_0
> 0$ in BE.

[Although $C(n)$ depends on $x_0$, we will see that the transport properties of the system are essentially independent of $x_0$. This explains why we do not write $C(n,x_0)$.]{.nodecor}

In terms of $C(n)$, one can write $$\begin{aligned}
L_{ij}^{(0)} &=& \frac{e^2}{h} C(0)  \Gamma_{ij}~,\label{L0}\\
L_{ij}^{(1)} &=&  \frac{e}{h} k_\B T C(1) \Gamma_{ij}~, \label{L1}\\
L_{ij}^{(2)} &=& \frac{(k_\B T)^2}{h} C(2) \Gamma_{ij}~. \label{L2}\end{aligned}$$ In Appendix A, we show that $C(0)$, $C(1)$ and $C(2)$ are positive for all $$x_0 = -\frac{\mu}{k_\B T} \in \left\{\begin{array}{l}(-1,\infty)
\mbox{ in MB}\\ 
(-\infty,\infty) \mbox{ in FD}\\ (0,\infty) \mbox{ in BE}\end{array}\right.~.$$ In the Maxwell-Boltzmann case, when $x_0 \in (-\infty,-1]$, we find that $C^{\rm MB}(0)$ and $C^{\rm MB}(2)$ are positive, but $C^{\rm
  MB}(1)$ is non-positive. As one may easily check, all the subsequent results also hold in this situation. Nevertheless, in order to be able to treat the three considered cases on the same footing we make the following restricting assumption.

In the Maxwell-Boltzmann case: $x_0 = -\mu/(k_\B T)\in (-1,\infty)$.\[Assumption xo\]

[Note that the coefficient $C^{\rm BE}(0)$ diverges as $x_0 \rightarrow 0$. This may be interpreted as a Bose-Einstein condensation, and shows that our theory breaks down in describing properly the transport of Bose-Einstein condensates.]{.nodecor}

\[remarks FD\] [In mesoscopic physics, one is usually interested in electronic transport at low temperature, $x_0
  = -\mu/(k_\B T) \rightarrow -\infty$. Observe that in this limit $C^{\rm FD}(1) \rightarrow 0$ and consequently $L_{ij}^{{\rm
FD}(1)} \rightarrow 0$. In order to have non-zero transport coefficients $L_{ij}^{{\rm
  FD}(1)}$ in this regime, one usually considers the first order term in the Taylor expansion, $\Gamma_{ij}(E) = \Gamma_{ij}(\mu) +
  \Gamma'_{ij}(\mu) (E-\mu) + O((E-\mu)^2)$, which leads to $L_{ij}^{{\rm
FD}(1)} \sim \Gamma'_{ij}(\mu)$.]{.nodecor}

[Observe that the transport coefficients $L_{ij}^{(k)}$ depend on the properties of the quantum dots through the scattering matrix $S$. Here we can see why the assumption A\[Assumption Gamma2\] is useful. Indeed, suppose that $t_{ii} = M_i$ for some $i$ (total reflection). Then $t_{ij} = 0$ for all $j \not = i$ and consequently $L^{(k)}_{ij} = 0$ for all $j$. In other words, in such a situation the $i$-th reservoir is disconnected from the rest of the system.]{.nodecor}

Looking at the relations –, one sees that the transport coefficients are related as follows: $$L_{ij}^{(0)} = \frac{e^2}{h} C(0) \Gamma_{ij}~,
\hspace{5mm}  L_{ij}^{(1)} = Q_1 L_{ij}^{(0)}  \hspace{5mm}
\mbox{and} \hspace{5mm}
L_{ij}^{(2)} = Q_2 L_{ij}^{(0)}~,$$ where $$Q_1 =\frac{k_\B T}{e} \frac{C(1)}{C(0)} > 0 \hspace{5mm}
\mbox{and} \hspace{5mm} Q_2 = \frac{k_\B^2 T^2}{e^2} \frac{C(2)}{C(0)} > 0~.\label{L4L1}$$

\[Autre remarks FD2\] [The quantity $Q_1$ is related to the Seebeck coefficient (see  \[Rem Seebeck\]). In the Fermi-Dirac case, $\lim_{x_0 \rightarrow  -\infty} Q_2/T^2 = (\pi^2/3) (k_\B/e)^2$ is called the Lorentz number.]{.nodecor}

Recalling that $\Gamma_{ij} =
M_i \delta_{ij} - t_{ij}$, one immediately obtains the following properties ($k = 0,1,2$): $$\label{Property Lij}
L^{(k)}_{ij} \ \left\{\begin{array}{ll} > 0 &
\mbox{if } \ i
= j\\
< 0 &
\mbox{if } \ i
\not = j
\end{array}\right.~,$$ and $$\label{Property Lij2}
\sum_{i} L^{(k)}_{ij} = 0~, \hspace{2mm} \forall j \hspace{7mm} \mbox{and} \hspace{7mm} \sum_{j} L^{(k)}_{ij} =
0~, \hspace{2mm} \forall i~.$$

As we shall see later on, the following ratio, $\mathcal{R}$, will play an important role. Using standard techniques, we show in Appendix B that (for all $x_0$):

$$\label{Ratio R}
\mathcal{R} \equiv \frac{Q_1}{\sqrt{Q_2}} =
\frac{C(1)}{\sqrt{C(0) \cdot C(2)}} \in (0,1)~.$$

[In terms of the matrix $$C = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}C(0) & C(1)\\C(1) & C(2)\\ \end{array}\hspace{-2mm}\right)},$$ the inequality $0 < \mathcal{R} < 1$ is equivalent to $\det(C) = C(0)
C(2) - C(1)^2 > 0$.]{.nodecor}

The Entropy Production {#The entropy production}
----------------------

In this subsection, we show that the transport matrix $L$ is real positive semi-definite. In other words, we show that the entropy production rate is non-negative: $$\label{entr prod}
\sigma_s = \sum_{i,j=1}^{2N+4}  L_{ij} V_i V_j \geq 0~,$$ where the thermodynamic forces are $$V_i = \left\{ \begin{array}{ll}\delta \mu_{i-1}/e & \mbox{ if } \; \, i = 1,
  \dots, N+2\\ \delta T_{i-(N+3)}/T &  \mbox{ if } \; \; i = N+3, \dots, 2N+4\end{array}\right. ~.$$ Here, we have set $\L = 0$ and $\R = N+1$. Let $X_i = \delta \mu_i/e$ and $Z_i = \sqrt{Q_2} \ \delta
T_i/T$, for $i \in \{0, 1, \dots, N, N+1\}$. Then, in Appendix C, we show that one can rewrite the entropy production rate as follows: $$\label{sigmas expression}
\sigma_s = \sum_{\scriptsize\begin{array}{c}
    i,j = 0\\ i < j \end{array}}^{N+1} (-L^{(0)}_{ij}) \ I_{ij}~,$$ where $$I_{ij} = (X_i - X_j)^2
  + (Z_i - Z_j)^2 -  2 \mathcal{R} C_{ij}~, \hspace{5mm} C_{ij} = X_i Z_j + X_j Z_i - X_i Z_i - X_j Z_j~.$$ Observe now that $(-L^{(0)}_{ij}) > 0$ for all $i \not = j$. Therefore, to show that $\sigma_s \geq 0$ it is sufficient to show that $I_{ij} \geq 0$ for all $i \not = j$. As we shall see, it will be crucial that the ratio $\mathcal{R}$ satisfies $0 < \mathcal{R} < 1$. Assume first that $C_{ij} \leq 0$, then ($
\mathcal{R} > 0$) $$\label{In1a}
I_{ij} \geq (X_i - X_j)^2 + (Z_i - Z_j)^2 \geq 0~.$$ Assume next that $C_{ij} > 0$, then ($\mathcal{R} < 1$) $$I_{ij} >  (X_i - X_j)^2 + (Z_i - Z_j)^2  - 2 C_{ij} = (X_i - X_j + Z_i - Z_j)^2 \geq 0~.$$ This ends the proof that $\sigma_s \geq 0$.

[The term $(X_i - X_j)^2$ accounts for the entropy production due to the electric current, and it is well known in mesoscopic physics (see [@Buttiker-Symmetry]). The term $(Z_i - Z_j)^2$ is related to the heat current, and one sees that there is also a thermoelectric term, $-2 \mathcal{R} C_{ij}$, which might take positive and negative values.]{.nodecor}

Equilibrium and Non-Equilibrium States {#Equilibrium
  and non-equilibrium states}
--------------------------------------

Here we present some equivalent characterizations of the equilibrium and non-equilibrium states of the multi-terminal system. We start with the following.

We say that the system is at *equilibrium* if $$T_\L = T_1 = \dots = T_N = T_\R  \hspace{5mm}\mbox{ and } \hspace{5mm}
\mu_\L = \mu_1 = \dots = \mu_N = \mu_\R~.$$ Otherwise, we say that the system is *out of equilibrium*.

Now, as one expects, the system is at equilibrium only in the situations in which all the electric and heat currents vanish, $$\{\mbox{System is at equilibrium}\} \iff \{I_i = 0 \hspace{3mm}
\mbox{and} \hspace{3mm} J_i = 0, \hspace{3mm} \forall i\}~.$$ This equivalence is obtained by writing the no-current condition, $I_i
= J_i = 0$ for all $i \in \{\L, \R, 1,\dots,N\}$, in a matrix form and by using the properties of the transport coefficients $L^{(k)}_{ij}$ as well as the fact that $0 < \mathcal{R} < 1$. The details are presented in Appendix D.

Also, one can characterize the equilibrium and non-equilibrium states of the system in terms of the entropy production: $$\left\{\hspace{-1mm}\begin{array}{cc}\mbox{System is at}\\\mbox{equilibrium}\end{array}\hspace{-1mm}\right\} \iff \sigma_s = 0 \hspace{5mm} \mbox{and}
\hspace{5mm} \left\{\hspace{-1mm}\begin{array}{cc}\mbox{System is}\\\mbox{out of
  equilibrium} \end{array}\hspace{-1mm}\right\}\iff \sigma_s > 0~.$$ These equivalences easily follow from the expression for $\sigma_s$ and the inequality . The details can be found in Appendix D.

The Self-Consistency Condition {#The null currents condition}
==============================

We now turn to the resolution of the self-consistency condition . Using the expressions – one can rewrite the self-consistency condition as follows ($i = 1,\dots,N$): $$\begin{aligned}
\sum_{j = 1}^{N} \left(L_{ij}^{(0)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(1)}
\, \frac{\delta T_j}{T}\right) &=& - \sum_{j = \L,\R} \left(L_{ij}^{(0)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(1)}
\, \frac{\delta T_j}{T}\right)~,\label{Condition1}\\
\sum_{j = 1}^{N} \left(L_{ij}^{(1)} \, \frac{\delta \mu_j}{e} + 
L_{ij}^{(2)}
\, \frac{\delta T_j}{T}\right) &=& - \sum_{j = \L,\R} \left(L_{ij}^{(1)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(2)}
\, \frac{\delta T_j}{T}\right)~.\label{Condition2}\end{aligned}$$ We recall that we are given the values $(T_\L, \mu_\L)$ and $(T_\R, \mu_\R)$, so that the right-hand side of the above equations are supposed to be known ($T_j = T + \delta T_j$ and $\mu_j = \mu + \delta \mu_j$). Therefore, the self-consistency condition constitutes a set of $2N$ equations for the $2N$ unknown variables $T_1, \dots, T_N$ and $\mu_1, \dots,
\mu_N$.

To solve these equations, it is convenient to introduce the vectors $X, Y \in \mathbb{R}^{N}$ and the $N\times N$ matrix $L^{(0)}_{\rm C}$ defined by ($i,j = 1,\dots,N$): $$X_j = \frac{\delta \mu_j}{e}~, \hspace{5mm} Y_j = \frac{\delta
  T_j}{T}~,  \hspace{5mm}
  \mbox{and} \hspace{5mm} (L^{(0)}_{\rm C})_{ij} = L^{(0)}_{ij}~.$$

[Recalling that $L^{(0)}$ is the $(N+2) \times (N+2)$ matrix defined in  \[def transport matrix\], one sees that $L^{(0)}_{\rm
  C}$ is the reduced matrix associated to the “central” terminals $1,\dots,N$.]{.nodecor}

We also define for $\ell = \L, \R$: $$X_\ell = \frac{\delta \mu_\ell}{e}~, \hspace{5mm} Y_\ell = \frac{\delta
  T_\ell}{T}~, \hspace{5mm}
  \mbox{and} \hspace{5mm} D_\ell = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}L^{(0)}_{1 \ell} \\ \vdots
  \\ L^{(0)}_{N \ell}\\ \end{array}\hspace{-2mm}\right)}~.$$ Then the equations – can be rewritten as follows: $$\begin{aligned}
L^{(0)}_{\rm C} \left(X + Q_1 Y\right) &=& - \sum_{\ell = \L, \R}
\left(X_\ell +  Q_1 Y_\ell \right) \ D_\ell~, \label{Eq XY1}\\
L^{(0)}_{\rm C} \left(Q_1 X + Q_2 Y\right) &=& - \sum_{\ell = \L, \R}
\left(Q_1 X_\ell +  Q_2 Y_\ell \right) \ D_\ell~.\label{Eq XY2}\end{aligned}$$ Since $\mathcal{R} \not = 1$, one has $Q_2 \not =
(Q_1)^2$. As a consequence, if one takes $(Q_2/Q_1)$ times and subtract $\eref{Eq XY2}$ and similarly if one takes $(-1/Q_1)$ times and add , one obtains $$L^{(0)}_{\rm C} X = - \sum_{\ell = \L, \R} D_\ell X_\ell \hspace{5mm}
  \mbox{and} \hspace{5mm} L^{(0)}_{\rm C} Y = - \sum_{\ell = \L, \R} D_\ell
  Y_\ell~. \label{Equ MXY}$$ We see here that the equations – decouple and that the equations for the chemical potentials and for the temperatures are actually identical. In order to proceed, we need to know that the matrix $L^{(0)}_{\rm C}$ is invertible. This fact is proved in Appendix E by showing that $L^{(0)}_{\rm C}$ is real positive definite. In particular, this shows that the self-consistency condition has a unique solution.

By inverting the matrix $L^{(0)}_{\rm C}$, we obtain the following self-consistent temperatures and chemical potentials ($i = 1,\dots,N$): $$\begin{aligned}
T_i &=& T_\L +  A_i \,
(T_\R - T_\L)~,\label{Eq2 Profile}\\
\mu_i &=& \mu_\L +  A_i \,
(\mu_\R - \mu_\L)~,\label{Eq1 Profile}\end{aligned}$$ where $$A_i = \sum_{j = 1}^{N} (\Gamma_{\rm C}^{-1})_{ij} \ t_{j \R}~, \hspace{5mm}
\mbox{with}  \hspace{5mm} \Gamma_{\rm C} = (\Gamma_{ij})_{i,j = 1}^{N}~. \label{Eq2 Profileb}$$ The details about their derivation are given in Appendix F.

By substituting these expressions into the relations – for the currents, with $i = \L$, and by using the properties of the transport coefficients $L^{(k)}_{ij}$, one deduces the (average) currents flowing through the system (from left to right): $$\begin{aligned}
I_\L &=& \sigma_0 \left( \frac{\mu_\R - \mu_\L}{e}\right) + \sigma_1
\left( \frac{T_\R - T_\L}{T}\right)~,\label{Current I}\\
J_\L &=& \sigma_1 \left(\frac{\mu_\R - \mu_\L}{e}\right) + \sigma_2  \left(\frac{T_\R - T_\L}{T}\right)~,\label{Current J}\end{aligned}$$ where $$\sigma_0 = L^{(0)}_{\L\R} + \sum_{j = 1}^{N} A_j \ L^{(0)}_{\L j}~,
\hspace{5mm} \sigma_1 = Q_1 \sigma_0 \hspace{5mm} \mbox{and} \hspace{5mm} \sigma_2 = Q_2 \sigma_0~. \label{Sigma1}$$

[The corresponding expressions to – in the case of electronic transport at low temperature ($T_i = 0$ and $J_i = 0$ for $i \in \{\L, \R, 1, \dots N\}$) were obtained in [@Amato].]{.nodecor}

\[Rem Seebeck\] [Setting $I_\L = 0$, one deduces the thermoelectric field $\mathcal{E} = - \nabla V =
  \mathcal{S} \nabla T$, where $\mathcal{S} = Q_1/T = k_\B/e \cdot
  C(1)/C(0) > 0$ is the thermopower or Seebeck coefficient. Note that $\mathcal{S}$ does not depend on $N$. For further discussions see [@Streda; @BeenakkerThermoPower; @StaringOscillations; @Molenkamp; @Godijn; @Mathias; @Kato].]{.nodecor}

Using the relations –, one may rewrite the currents – as follows: $$\begin{aligned}
I_\L &=& \sum_{j} L^{(0)}_{\L j}  \left( \frac{\mu_j -
  \mu_\L}{e}\right) +  L^{(1)}_{\L j} \left(\frac{T_j - T_\L}{T}\right)~,\\
J_\L &=& \sum_{j} L^{(1)}_{\L j}  \left( \frac{\mu_j -
  \mu_\L}{e}\right) +  L^{(2)}_{\L j} \left(\frac{T_j - T_\L}{T}\right)~.\end{aligned}$$

Although the multi-terminal system considered so far has no specific geometry (see  \[The General System\]), we shall nevertheless use the word “profile” for the arrangements $\mu_1, \dots,
  \mu_N$ and $T_1, \dots, T_N$. This terminology will be fully justified in  \[Numerical analysis\], where we will restrict our attention to the linear chain of quantum dots ( \[The Quantum System\]).

[Note that the temperature and chemical potential profiles are decoupled, $T_i$ only depends on $(T_\L, T_\R)$ and $\mu_i$ only depends on $(\mu_\L, \mu_\R)$. On the other hand, since $Q_1 > 0$, there can be a heat current due to a gradient of chemical potential ($T_\L = T_\R$), and reciprocally, there can be an electric current due to a gradient of temperature ($\mu_\L =
  \mu_\R$).]{.nodecor}

\[remark on B\] [We see that the form of the profiles is given by the coefficients $A_1, \dots, A_N$. Hence, in what follows we will use the word “profile” to refer to these coefficients. Looking at the expression for $A_i$ one sees that it does not depend on the distribution function – describing the reservoirs. In this sense, we say that the profile is *universal*. One may wonder whether this is true for any distribution function $f$. It seems not to be the case. Indeed, the key ingredients to obtain the preceding results are: (i) the specific form of $f$ and (ii) the ratio $\mathcal{R}$, defined in , must satisfy $0 < \mathcal{R} <
1$. Basically, we expect that any distribution function $f$ satisfying (i)-(ii) will lead to the same results. Note that universal temperature profiles were also found for a quantum harmonic chain coupled at both ends to two phonon reservoirs at different temperatures [@Saito].]{.nodecor}

\[remark on B\] [Although the profiles are the same in the three considered situations, one can distinguish between the classical and quantum chain of quantum dots. Indeed, assume we are given the scattering matrices, $S^{(1)}, \dots, S^{(N)}$, associated to the $N$ quantum dots. Then, as explained in detail in Appendix H, one can compose them into a global scattering matrix $S$, out of which one extracts the transmission probabilities $t_{ij}$ (see ). This is the natural *quantum* way of working.]{.nodecor}

[In order to extract the interference effects, one may first compute the probability matrices, $P^{(1)}, \dots, P^{(N)}$, defined by $P^{(k)}_{ij;mn} = |S^{(k)}_{ij;mn}|^2$, then compose them into a global probability matrix $P$ and finally set $t_{ij} = \sum_{m=1}^{M_i}
\sum_{n=1}^{M_j} P_{ij;mn}$. This way of computing $t_{ij}$ will be referred to as *classical*.]{.nodecor}

[As one may check, the composition law of the probability matrices is the same as for the scattering matrices, all expressions in Appendix H hold with $P$ instead of $S$. In this sense, we will study in  \[Numerical analysis\] the classical versus quantum situations (see also [@ShapiroB; @CahayArray]).]{.nodecor}

In Appendix G, we derive an interesting alternative expression for the coefficients $A_i$. We obtain ($i = 1, \dots, N$) $$\label{Prop B}
A_i = \frac{\sum_{j = 1}^{N} (-1)^{i+j} \det\left(\Gamma_{\rm C}(j,i)\right)
  \ t_{j \R}}{\sum_{j = 1}^{N} (-1)^{i+j} \det\left(\Gamma_{\rm C}(j,i)\right)
  \ [t_{j \L} + t_{j \R}]}~,$$ where $\Gamma_{\rm C}(j,i)$ denotes the $(j,i)$ minor of $\Gamma_{\rm  C}$.

\[Remark on B1\] [As an illustration, let us apply the relation in the case $N = 1$. One has $$\label{Eq on B1}
A_1 = \frac{t_{1 \R}}{t_{1 \L} + t_{1 \R}} \hspace{5mm}
\Longrightarrow \hspace{5mm}\mu_1 = \frac{t_{1\L} \ \mu_\L + t_{1\R} \
  \mu_\R}{t_{1 \L} + t_{1 \R}}~.$$ This relation for $\mu_1$ corresponds to the one obtained in [@Buttiker-Symmetry] for a three-terminal conductor under the conditions $T_\L = T_1 = T_\R =
0$ (which implies $J_1 = 0$, since $\delta T_j = 0$ and $L_{ij}^{{\rm FD} (1)} = 0$ in this case) and $I_1 = 0$.]{.nodecor}

An important consequence of is ($i = 1, \dots, N$): $$\label{Proposition B Sigma}
A_i \in (0,1)~.$$ This implies that the temperatures and chemical potentials of the reservoirs, connected self-consistently to the system, are bounded between the boundary values $T_\L, T_\R$ and $\mu_\L, \mu_\R$, respectively. In other words, for $i =
  1, \dots,N$, one has $$\min\{ T_\L, T_\R \} \leq T_i \leq \max\{ T_\L, T_\R \}
\hspace{5mm}\mbox{and} \hspace{5mm} \min\{ \mu_\L, \mu_\R \} \leq \mu_i \leq \max\{ \mu_\L, \mu_\R \}~.$$

[The property also implies that the self-consistent parameters $T_1, \dots, T_N$ and $\mu_1, \dots, \mu_N$ are increasing functions of $T_\L$, $T_\R$, and $\mu_\L$, $\mu_\R$, respectively.]{.nodecor}

Moreover, looking at the expressions for the coefficients $\sigma_0$, $\sigma_1$ and $\sigma_2$, and recalling that $L^{(0)}_{ij} < 0$ if $i \not = j$, $Q_1 > 0$ and $Q_2 > 0$, one deduces that for $k = 0,1,2$: $$\label{Proposition B Sigma2}
\sigma_k < 0~.$$

[There are some simple consequences of with respect to –. (i) If $\mu_\L = \mu_\R$, then the heat current flows from the hot reservoir to the cold one, and similarly (ii) if $T_\L = T_\R$, then the electric current goes from the high chemical potential to the lower one. Note, however, that by setting $\mu_\L, \mu_\R$ and $T_\L, T_\R$ to appropriate values there might be a heat current going from the cold reservoir to the hot one. (A similar statement holds for the electric current.)]{.nodecor}

Ohm and Fourier Laws {#Ohm and Fourier laws theory}
--------------------

From the relations – for the electric and heat currents, one can find explicit expressions for the electric and heat conductances. If $T_\L = T_\R$, then setting $\mu = e V$, where $V$ denotes the electric potential, one obtains $$\label{Eq Ohm global1}
I_\L = - G_{\rm e} \ (V_\R - V_\L)~,$$ where the electric conductance $G_{\rm e}$ is given by $$G_{\rm e} = - \sigma_0~.$$ If $I_\L = 0$, one deduces $$\label{Eq Fourier global1}
J_\L = - G_{\rm h} \ (T_\R - T_\L)~,$$ where the heat conductance $G_{\rm h}$ is given by $$G_{\rm h} = \frac{\sigma_1^2 - \sigma_0 \sigma_2}{\sigma_0 T} = -
\frac{(Q_2 - Q_1^2)}{T} \ \sigma_0 = - \frac{k_\B^2 T}{e^2} \ \frac{C(0)
  C(2) - C(1)^2}{C(0)^2} \ \sigma_0~.$$ Recalling that $0 < \mathcal{R} < 1$, one sees that $Q_2 >
(Q_1)^2$. Consequently, since $\sigma_0 < 0$, one deduces that the electric and heat conductances are positive: $$G_{\rm e} > 0 \hspace{5mm}  \mbox{and} \hspace{5mm} 
G_{\rm h} > 0~.$$

Assume now that the system has the linear geometry represented in  \[The Quantum System\] and let $a$ denote the length of the quantum dots and $\Sigma$ denote the width of the openings connecting any two dots (the effective cross section of the system). Since $L=Na$ corresponds to the length of the system made of $N$ quantum dots, one may introduce the global gradients: $$\label{global gradients}
\nabla V = \frac{V_\R
- V_\L}{L} \hspace{7mm} \mbox{and}  \hspace{7mm} \nabla T = \frac{T_\R
- T_\L}{L}~.$$ Then, under the same conditions as in and , the current *densities*, $\mathcal{I}_\L = I_\L /\Sigma$ and $\mathcal{J}_\L = J_\L /\Sigma$, are given by the *global* Ohm and Fourier laws: $$\label{Eqs global}
\mathcal{I}_\L = - \kappa_{\rm e} \nabla V \hspace{5mm} \mbox{and} \hspace{5mm} \mathcal{J}_\L = - \kappa_{\rm h} \nabla T~,$$ where the electric and heat conductivities, $\kappa_{\rm e}$ and $\kappa_{\rm h}$, are given by $$\label{Eqs cond global}
\kappa_{\rm e} = \frac{L G_{\rm e}}{\Sigma} > 0 \hspace{5mm} \mbox{and} \hspace{5mm}
\kappa_{\rm h} = \frac{L G_{\rm h}}{\Sigma} > 0~.$$

[In the Fermi-Dirac situation, one has $$\lim_{x_0 \rightarrow 0} Q_1 = 0 \hspace{2mm} \mbox{and}
  \hspace{2mm} \lim_{x_0
  \rightarrow 0} Q_2 = \frac{\pi^2}{3} \left(\frac{k_\B}{e}\right)^2 T^2
  \hspace{2mm} \Longrightarrow  \hspace{2mm} \lim_{x_0 \rightarrow 0} \frac{\kappa_{\rm h}}{\kappa_{\rm e} T} = \frac{\pi^2}{3} \left(\frac{k_\B}{e}\right)^2~.$$ The last relation is the Wiedemann-Franz law giving the Lorentz number.]{.nodecor}

In a similar manner, one may introduce the following local gradients:

$$\nabla V(i) = \frac{V_{i+1} - V_i}{a} \hspace{5mm} \mbox{and} \hspace{5mm}
\nabla T(i) = \frac{T_{i+1} - T_i}{a}~.$$

Let $\mathcal{I}(i) = \mathcal{I}_\L$ and $\mathcal{J}(i) =
 \mathcal{J}_\L$, where $i \in \{1, \dots, N-1\}$, denote the net local current densities flowing from the $i$-th quantum dot into the $(i+1)$-th quantum dot. Then, using the equations – and , one obtains the *local* Ohm and Fourier laws: $$\label{Local OF}
\mathcal{I}(i) = - \kappa_{\rm e}(i) \nabla V(i) \hspace{5mm} \mbox{and}
\hspace{5mm} \mathcal{J}(i) = - \kappa_{\rm h}(i) \nabla T(i)~,$$ where the *local* electric and heat conductivities are given by $$\kappa_{\rm e}(i) = \frac{\kappa_{\rm e}}{N(A_{i+1} - A_i)} \hspace{5mm}
\mbox{and} \hspace{5mm} \kappa_{\rm h}(i) = \frac{\kappa_{\rm h}}{N(A_{i+1} - A_i)}~.$$

[Although the global conductivities, $\kappa_{\rm e}$ and $\kappa_{\rm h}$, are always positive, we will see, in  \[Numerical analysis\], that the local conductivities, $\kappa_{\rm e}(i)$ and $\kappa_{\rm h}(i)$, may take negative values.]{.nodecor}

The coefficients $\sigma_0$, $\sigma_1$ and $\sigma_2$, the conductances $G_{\rm e}$ and $G_{\rm h}$, and the conductivities $\kappa_{\rm e}$ and $\kappa_{\rm h}$ depend on the distribution function $f$ describing the reservoirs. Nevertheless, they are all proportional to $\sigma_0$, with some multiplicative factor independent of $N$. The important observation is that the dependence on the nature of the particles is entirely contained in these multiplicative factors. This permits to introduce the following universal quantity.

\[universal conductance\] We define the *universal conductivity* as $$\label{universal conductance1}
\kappa(N) = -\frac{h}{e^2 C(0)} N \sigma_0  = N \left[t_{\L\R} +
  \sum_{i = 1}^{N} \sum_{j = 1}^{N} t_{\L j} \ (\Gamma_{\rm
  C}^{-1})_{ji} \ t_{i \R}\right]~.$$

\[Remark on Sigma1\] [As an illustration, let us consider the case $N=1$. One has $$\label{equ kappa1}
\kappa(1) = t_{\L \R} + \frac{t_{\L 1} \ t_{1 \R}}{t_{1 \L} + t_{1 \R}}~.$$ As expected, the universal conductivity is the sum of two contributions: a direct coherent transmission and an indirect incoherent transmission through the intermediate terminal. This result is well known in mesoscopic physics [@Buttiker-Symmetry], and one sees that the relation is the generalisation to $N$ self-consistent reservoirs.]{.nodecor}

Looking at the expression for the universal conductivity, one sees that $\kappa(N) > 0$, since $\sigma_0 <
0$ and $C(0) > 0$. Using moreover the relations and $t_{\L \R} < \min\{M_\L,
M_\R\}$, one finds the following *optimal* bounds: $$0 < \kappa(N) < N \cdot \min\{M_\L, M_\R\}~.\label{bounds}$$

Numerical Simulations {#Numerical analysis}
=====================

The Chain of Quantum Dots
-------------------------

To prove the validity of the *global* Ohm and Fourier laws in our model, it remains to show that the transport is *normal*, that is, when the system becomes macroscopic the global electric and heat conductivities, $\kappa_{\rm e}$ and $\kappa_{\rm h}$, become intrinsic properties of the macroscopic system and consequently become independent of its length $L = N a$. More precisely, one has to prove that the global electric and heat conductivities converge to finite limits as $L \rightarrow \infty$. Looking at the relations and , one sees that $\kappa_{\rm e}$ and $\kappa_{\rm h}$ are both proportional to $\kappa(N)$, with some multiplicative factor independent of $N$, and consequently it is sufficient to investigate the limit of $\kappa(N)$ as $N \rightarrow \infty$.

As explained previously, the multi-terminal system considered so far has no specific geometry ( \[The General System\]), and therefore $\kappa(N)$ will not converge as $N \rightarrow~\infty$ in general. In order obtain a finite limit for $\kappa(N)$ as $N \rightarrow
\infty$ we shall proceed as follows.

Let $\left(S^{(k)}\right)_{k = 1}^{\infty}$ be an infinite sequence of scattering matrices, each one being associated to a quantum dot, and let $S_N$ denote the *composite* global scattering matrix associated to $\left(S^{(1)},\dots,S^{(N)}\right)$. The details concerning the construction of $S_N$ are given in Appendix H. Having set-up the linear geometry represented in  \[The Quantum System\], one may expect that $\kappa(N)$ will converge to a finite value, $\kappa^{\infty}$, as $N
\rightarrow \infty$, and this is what we found in *typical* numerical experiments (see below). Unfortunately, the composite scattering matrix, $S_N$, is given in a rather cumbersome form in terms of $S^{(1)},\dots,S^{(N)}$, which prevents us from proving the existence of $\lim_{N \rightarrow \infty}
\kappa(N)$.

Nevertheless, we will see by using numerical simulations under random matrix theory, that the universal conductivity $\kappa(N)$ typically admits a finite limit, $\langle  \kappa^{\infty}  \rangle$, as $N \rightarrow
\infty$. This shows that the *global* Ohm and Fourier laws hold on statistical average, and we may write $$\label{eq Ge}
\langle \kappa^{\infty}_{\rm e}  \rangle = \frac{a}{\Sigma} \ \frac{e^2 C(0)}{h} \ \langle \kappa^{\infty}  \rangle$$ and $$\label{eq Gh}
\langle  \kappa^{\infty}_{\rm h}  \rangle =  \frac{a}{\Sigma} \ \frac{k_\B^2
  T}{h} \ \frac{C(0) C(2) - C(1)^2}{C(0)} \
\langle  \kappa^{\infty}  \rangle~.$$

Moreover, we will see that typically the profile, $A_1,\dots,A_N$, becomes linear as the number of dots $N \rightarrow
\infty$. Introducing the variable $x = i/(N+1) \in [0,1]$, this means that in the thermodynamic limit $N \rightarrow
\infty$, one has $$\langle 
A(x)  \rangle = x~,$$ and consequently $$\langle V(x) \rangle = V_\L + (V_\R - V_\L) \ x \hspace{5mm} \mbox{and}  \hspace{5mm}
\langle T(x) \rangle = T_\L + (T_\R - T_\L) \ x~.$$ Hence, the *local* Ohm and Fourier laws also hold on statistical average, with position-independent electric and heat conductivities.

[Here we compare the values of the average electric and heat conductivities for the three considered situations. Indeed, looking at the expressions –, one sees that they only differ from a multiplicative factor, which depends on the nature of the particles. In the electric and heat cases, we found the following order: BE $>$ MB $>$ FD, for all $x_0 > 0$ (the common admissible region).]{.nodecor}

Random Matrix Theory
--------------------

As mentioned in the introduction, the classical EY-model describing particle and energy transport ( \[The Classical System\]) has the following important property: it is chaotic. More precisely, the cells, *without* the discs, which are making up the system are ergodic. In analogy, we assume the following.

The quantum dots are classically chaotic.

Under this assumption, it is natural to describe the transport properties of the quantum dots with random matrix theory (RMT) [@Beenakker]. More precisely, we will consider that the quantum dots are described by the Wigner-Dyson circular ensembles (with $\beta = 1, 2$).

Here, we restrict our attention to the one channel situation,  we assume that $M_i = 1$ for $i \in \{\L, \R, 1, \dots,
N\}$. Now we consider that the $3 \times 3$ complex random matrices, $S^{(1)},\dots,S^{(N)}$, are independent and identically distributed over $U(3)$ with some measure. More specifically, we consider the following well-known cases:

-   The Circular Orthogonal Ensemble (COE, $\beta = 1$) describes the case of time-reversible systems.

-   The Circular Unitary Ensemble (CUE, $\beta = 2$) corresponds to the case in which the time-reversal symmetry is broken. Such a situation appears, for example, when the system is submitted to some external magnetic field.

In the sequel, we shall be interested mainly in the following quantities:

-   The universal profile $A_1, \dots, A_N$, which gives the shape of the temperature and chemical potential profiles. In the subsequent plots, we naturally set $A_0 = A_\L = 0$ and $A_{N+1}=A_\R = 1$.

-   The universal conductivity $\kappa(N)$, which gives the electric and heat currents through the chain of $N$ quantum dots.

As explained previously, we shall compare the classical and quantum situations by working in terms of the probability matrices, $P^{(1)},\dots,P^{(N)}$, and scattering matrices, $S^{(1)},\dots,S^{(N)}$, respectively (see  \[remark on B\]).

For a review on RMT we refer to [@Beenakker]. To implement the numerical simulations, we have followed [@Mezzadri]. In the sequel, the averages are made over an ensemble of $10^4$-$10^6$ realizations, depending on the size $N$ of the system.

The Global Transmission Probabilities {#The S-matrix circular ensembles}
-------------------------------------

We denote by $\langle  t^{(k)}_{ij}  \rangle_{1}$ and $\langle 
t^{(k)}_{ij}  \rangle_{2}$ the average of the transmission probabilities in COE($\beta = 1$) and CUE ($\beta = 2$), respectively. It turns out that one can determine their exact values [@Beenakker]: $$\label{Eq Transmission RMT}
\langle t^{(k)}_{ij}  \rangle_{1} = \left\{\begin{array}{ll}
1/4 & \mbox{if} \hspace{3mm}  i \not = j\\
1/2 & \mbox{if} \hspace{3mm}  i =
j\end{array}\right. \hspace{5mm} \mbox{and} \hspace{5mm} \langle 
t^{(k)}_{ij}  \rangle_{2} = \frac{1}{3}~, \hspace{3mm} \forall i,j~.$$ Now we use the construction explained in Appendix H to build the scattering matrix $S$ associated to $N$ quantum dots. As the reader may guess, it is then very difficult to make any analytical statement concerning the averages $\langle t_{ij}  \rangle_{1}$ and $\langle  t_{ij}  \rangle_{2}$. We thus turned to numerical simulations, and made the following observations for large system sizes $N$:

-   The classical and quantum global transmission probabilities, $t^{\rm{cl}}_{ij}$ and $t^{\rm{qu}}_{ij}$, are *different*, but their averages, $\langle  t^{\rm{cl}}_{ij}
      \rangle_{\beta}$ and $\langle  t^{\rm{qu}}_{ij}
      \rangle_{\beta}$, are the same, showing that *on average* the interferences are negligible when $N$ is large.

-   The averages $\langle  t_{ij}  \rangle_{\beta}$ do not depend on the system size $N$.

-   They are symmetric: $\langle  t_{ij}  \rangle_\beta = \langle  t_{ji} 
    \rangle_\beta$.

-   The average couplings between the terminals, $\langle  t_{ij}  \rangle_{\beta}$, depend on $|i-j|$ and are *short ranged*, $\langle t_{i j} \rangle_{\beta} \simeq 0$ if $|i-j|
    > 2$. In particular, the average probability matrix has *essentially* the following form: $$\langle  P  \rangle_{\beta} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}* & * & * &  & & & &  \\  \hspace{0.001mm} * & * & * & &  & & &\\ &
      * & * & *  & & &  & \\ &  & * & * & * & & & \\ & & & \ddots & \ddots
      & \ddots & &
      \\ & & & & * & * & * &  \\  & & & & & * & * & *  \\
      & & & &  &  * & * & *\\ \end{array}\hspace{-2mm}\right)}~.
    \begin{picture}(0,0)
    \put(-62,37){\Huge{0}}
    \put(-128,-47){\Huge{0}}
    \end{picture}$$ The values of the matrix elements represented by 0 are smaller than 0.05.

A typical example illustrating the above observations is given in  \[Fig Smatrix\], where $N = 20$. In this figure, five terminals are considered: $\ell \in \{\L, 5, 10, 15, \R\}$ (with $\L = 0$ and $\R = 21$). The different curves correspond to $\langle  t_{i\ell} 
\rangle_{\beta}$ for some $i \in \{\ell - 4, \dots, \ell + 4\}$. For example, the middle curves give the following values: $\langle  t_{7, 10} 
\rangle_{\beta}$, $\dots$, $\langle  t_{10, 10} 
\rangle_{\beta}$, $\dots$, $\langle  t_{13, 10} 
\rangle_{\beta}$.

(0,0) (155,235)[COE]{} (385,235)[CUE]{}

The Profiles and Currents {#The profiles and conductances}
-------------------------

Let us recall the expressions for the coefficient $A_i$ and for the universal conductivity $\kappa(N)$: $$A_i = - \sum_{j = 1}^{N} (\Gamma_{\rm C}^{-1})_{ij} \ t_{j \R}
\hspace{5mm} \mbox{and} \hspace{5mm} \kappa(N) = N \left[t_{\L\R} +
  \sum_{i = 1}^{N} \sum_{j = 1}^{N} t_{\L j} \ (\Gamma_{\rm
  C}^{-1})_{ji} \ t_{i \R}\right]~.$$

In  \[Fig B profile\], we show the average profiles $\langle
A(x) \rangle_\beta$ for $N =
5$. If $T_\L < T_\R$, these curves correspond essentially to the temperature profiles along the system. (The same holds for the chemical potential profiles if $\mu_\L < \mu_\R$). We see that all these profiles are not linear, the non-linearity being more stressed in the quantum situation, showing the effect of the interferences. However, we found that all these profiles, while conserving their shape, get closer and closer to linear as the number of dots increases and that eventually they become linear in the limit $N \rightarrow \infty$.

(0,0) (30,205)[COE]{} (262,205)[CUE]{}

\[rem toy\] [In order to emphasise the role of interferences, let us consider the following simple toy models. We consider that the classical probability matrices are identical and given by $P^{(k)}_{ij} = 1/3$, for all $i,j$. In such a situation, the profile $A(x)$ is linear (for all $N$). On the other hand, the following unitary scattering matrices $$S^{(k)} = \frac{1}{\sqrt{3}} {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} 1 & e^{\frac{2\pi}{3}i} &
  e^{\frac{2\pi}{3}i}\\ e^{\frac{2\pi}{3}i} & 1 & e^{\frac{2\pi}{3}i}
  \\ e^{\frac{2\pi}{3}i} & e^{\frac{2\pi}{3}i} & 1 \\ \end{array}\hspace{-2mm}\right)}~,$$ which satisfy the equiprobability property, $|S^{(k)}_{ij}|^2 =
\frac{1}{3}$, do not lead to a linear profile, showing the interference effects (see  \[Fig B profile Special2\]). A discussion about the conductivities in these cases is given in  \[do not need chaos\].]{.nodecor}

\[remark toy model\] [One may have noticed that the average temperature profiles differ from those in the EY-model. Indeed, in the EY-model the temperature profiles were either linear, convex or concave functions of the cell position, depending on the values of the injection rates $\gamma_\L$ and $\gamma_R$. Basically, the reason was that, in the presence of a particle current, for example, going from right to left, the disc $i$ hears more often from cell $i+1$ than from cell $i-1$. It has therefore a greater tendency to equilibrate with the right than with the left, causing a curvature in the temperature profile. In particular, the temperature profile depends on all the boundary values $T_\L$, $T_\R$, $\gamma_\L$ and $\gamma_\R$. In the present model, the coupling among the terminals is governed by the scattering matrix $S$ and we have assumed it is energy-independent leading to a self-consistent temperature profile dependent only on $S$ and on the boundary values $T_\L$ and $T_\R$. However, we expect that in the general energy-dependent S-matrix situation, the temperature profile will also depend on $\mu_\L$ and $\mu_R$. To see whether in this case the EY-model and our model agree needs further investigations.]{.nodecor}

We have shown in the preceding section that, if $\mu_\L = \mu_\R$ and $T_\L < T_\R$, then the heat current $J_\L$ goes from the hot reservoir (Right) to the cold one (Left), as one expects from the second law of thermodynamics. In  \[Fig B profile\], one sees that the profiles are monotically increasing from left to right. This means that, on average, heat also flows locally from hot to cold, on statistical average, the second law of thermodynamics also holds locally.

All the preceding results concern statistical average values. In  \[Fig B profile Special\], we present different realizations and one sees that very different profiles may occur. Roughly speaking, by choosing an appropriate combination of the S-matrices, $S^{(1)}, \dots, S^{(N)}$, one may generate almost any profile.

In particular, one sees in  \[Fig B profile Special\] that there exist *particular* realizations for which the heat current flows *locally* from cold to hot. In other words, the profiles are not monotone increasing for every realization, but only on statistical average. Such phenomena are well known in mesoscopic systems (see  [@Buttiker-Symmetry; @ButtikerNR; @ButtikerPotential]).

We observed that the quantum dots about which such phenomena occur (such that $A_{i} < A_{i-1}$ and consequently $\kappa_{\rm e}(i-1) < 0$ and $\kappa_{\rm h}(i-1) < 0$) have typically a very high reflection probability $t_{ii}$. As the system size $N$ grows, we see that such phenomena diminish in intensity and frequency. Basically, we expect that, when the system is macroscopic, the second law of thermodynamics also holds locally in every realization.

Moreover, in the thermodynamic limit $N \rightarrow \infty$, we see that the average profiles become linear in all situations (COE, CUE, classical and quantum): $$\label{Profile Ax}
\langle A(x) \rangle_\beta = x~, \hspace{5mm} \forall x \in [0,1]~.$$ Consequently, $$\langle V(x) \rangle_\beta = V_\L + (V_\R - V_\L) \ x \hspace{5mm} \mbox{and}  \hspace{5mm}
\langle T(x) \rangle_\beta = T_\L + (T_\R - T_\L) \ x~.$$

We next present some results concerning the universal conductivity $\kappa(N)$. In  \[Fig Sigma Small N\], we show $\langle \kappa(N)
\rangle_\beta$ for small $N$. Observe that in the COE cases the conductivity decreases with $N$ while in the CUE cases it first increases and only then decreases.

For larger values of $N$, we see in  \[Fig Sigma\] that the average universal conductivity seems to admit a finite limit $\langle \kappa^{\infty} \rangle_\beta$ as $N
\rightarrow \infty$. To support this observation and obtain a numerical value for $\langle \kappa^{\infty} \rangle_\beta$, we have proceeded as follows. We found that the curves in  \[Fig Sigma\] are very well fitted (the largest chi-square being 0.035) by functions (not shown in the figures) of the following form: $\langle  \kappa^{\infty} 
\rangle_\beta + c_\beta/N^\alpha$. Here, $c_\beta$ is some positive constant and $\alpha$ is an exponent found to be the same in COE and CUE, and given by $\alpha^{\rm cl} = 1/2$ and $\alpha^{\rm qu} = 1$. Moreover, $$\begin{aligned}
\mbox{Classical:} & &  \hspace{10mm}\langle \kappa^{\infty} \rangle_1 = 0.057  \hspace{7mm} \mbox{and}
\hspace{7mm} \langle  \kappa^{\infty}  \rangle_2 = 0.531~.\\ 
\mbox{Quantum:} & & \hspace{10mm} \langle  \kappa^{\infty} \rangle_1 = 0.005  \hspace{7mm} \mbox{and}
\hspace{7mm} \langle  \kappa^{\infty}  \rangle_2 = 0.164~.\end{aligned}$$

(0,0) (158,195)[COE]{} (390,195)[CUE]{}

Let us denote by $\langle \kappa^{\infty}_{\rm e} \rangle_\beta$ and $\langle \kappa^{\infty}_{\rm h} \rangle_\beta$ the corresponding limits: $$\label{exp Ge}
\langle \kappa^{\infty}_{\rm e} \rangle_\beta = \frac{a}{\Sigma} \ \frac{e^2 C(0)}{h} \ \langle 
\kappa^{\infty}  \rangle_\beta$$ and $$\label{exp Gh}
\langle \kappa^{\infty}_{\rm h} \rangle_\beta  = \frac{a}{\Sigma} \ \frac{k_\B^2 T}{h} \ \frac{C(0) C(2) - C(1)^2}{C(0)} \
\langle  \kappa^{\infty}  \rangle_\beta~.$$ These numerical results show that the global Ohm and Fourier laws hold on *statistical average* in our chain of quantum dots with self-consistent reservoirs. Since the average temperature and chemical potential profiles become linear in the limit $N \rightarrow \infty$, it follows that the local Ohm and Fourier laws also hold on statistical average, with position-independent electric and heat conductivities, $$\langle \mathcal{I}(x) \rangle_\beta = - \langle \kappa^{\infty}_{\rm
  e}(x) \rangle _\beta \, \nabla \langle V(x) \rangle_\beta
\hspace{4mm}\mbox{and} \hspace{4mm} \langle \mathcal{J}(x)
\rangle_\beta = - \langle \kappa^{\infty}_{\rm h}(x) \rangle_\beta \, \nabla \langle T(x) \rangle_\beta ~,$$ with $\langle \kappa^{\infty}_{\rm e}(x)  \rangle_\beta = \langle \kappa^{\infty}_{\rm e} \rangle_\beta$ and $\langle \kappa^{\infty}_{\rm h}(x) \rangle_\beta = \langle \kappa^{\infty}_{\rm h} \rangle_\beta$ for all $x \in (0,1)$. Here we have kept the global gradients constant while taking the limit $N \rightarrow \infty$.

\[rem optimal\]

Although we have shown that typically (in the sense of RMT) $\kappa(N)$ converges to a finite limit as $N \rightarrow
  \infty$, one deduces from the examples in  \[rem optimal\] that the linear geometry, represented in  \[The Quantum System\], is not sufficient to ensure the existence of such a limit. The characterization of the sequences of scattering matrices, $\left(S^{(k)}\right)_{k = 1}^{\infty}$, for which the limit, $\lim_{N \rightarrow \infty} \kappa(N)$, exists and is finite remains an open problem. Note however that, approximate results were obtained in some particular cases [@Amato; @roy]. Note also that the particular ordered cases, in which all the local scattering matrices are identical, are not simpler to handle with the relations given in Appendix H.

In  \[Fig Sigma\], one sees two manifestations of weak localization:

-   The classical conductivity is larger than the quantum conductivity. Hence, the effect of the interferences is to decrease the intensity of the electric and heat currents. Note that this effect is rather subtle. Indeed, one has $$\langle \kappa(N)  \rangle_\beta = N \left[\langle  t_{\L\R}  \rangle_\beta
      +
      \sum_{i,j = 1}^{N} \langle  t_{\L j} \ 
      (\Gamma_{\rm C}^{-1})_{ji} \ t_{i \R}   \rangle_\beta\right]~.$$ Since the average global transmission probabilities are the same in the classical and quantum situations, the effect of the interferences on the conductivity must be due to correlations among the random variables $t_{ij}$.

-   The application of an external field (CUE) increases the values of the conductivities and consequently increases the intensity of the electric and heat currents. This can be understood as follows: In equation  and in  \[Fig Smatrix\], one sees that the presence of a magnetic field will (on average) increase the transmissions probabilities.

\[Remark Amato\] [The case of electronic transport at low temperature was already considered by D’Amato and Pastawski [@Amato]. Modelling the system as a nearest-neighbour tight-binding Hamiltonian they found *approximate* expressions for the global transmission probabilities $t_{ij}$, which turn out to satisfy the properties (2)–(4) presented in  \[The S-matrix circular ensembles\]. From these approximate $t_{ij}$ they deduced that, in the limit the system size $N \rightarrow \infty$, the self-consistent chemical potential profile is linear and the universal conductivity is finite, in agreement with our results. However, in our opinion, these (nice) results need further analysis to constitute a rigorous derivation of Ohm’s law in their model.]{.nodecor}

Chaos, Disorder and Decoherence {#do not need chaos}
-------------------------------

In this paper, we were not interested in finding the optimal conditions to obtain Ohm and Fourier laws, but rather to establish a quantum version of the classical chaotic EY-model. This led us naturally to use RMT to characterize the quantum dots. Nevertheless, one may wonder whether chaoticity or randomness, as modelled by RMT, is essential to obtain normal transport, for Ohm and Fourier laws to hold. The answer is no, since, for example, the two toy models presented in  \[rem toy\], which may be interpreted as ordered non-chaotic cases, do lead to Ohm and Fourier laws. Indeed, one sees in  \[Conductivity toy model\] that their corresponding universal conductivities converge to finite values as the number of dots increases. Observe however that, contrary to the RMT situations, the conductivities increase with the number of dots and the “quantum” universal conductivity is higher than the “classical” one. This shows in particular that, although interferences are typically (in the sense of RMT) destructive, they may be constructive in some particular cases. Actually, this phenomenon is typical in *ordered* chains (see below).

Another natural investigation concerns the effects of disorder. We say that the chain is *ordered* if all the local scattering matrices are identical ($S^{(1)} = \dots = S^{(N)}$) and *disordered* otherwise. The results presented in  \[Fig B profile\] and \[Fig Sigma\] concern the disordered chaotic cases and in  \[Fig B profile Special2\] and \[Conductivity toy model\] one can find two examples of ordered non-chaotic cases. The ordered chaotic situation is obtained by using RMT with the condition that for each realization $S^{(1)} = \dots = S^{(N)}$. The results are given in  \[Fig B profile ordered\] and \[Conductivity Ordered model\]. Observe that by applying a magnetic field the “s-shapes” of the profiles are reversed and that the conductivities behave qualitatively as in the (ordered) toy model cases. Note also that in the quantum COE situation, the asymptotic value of the universal conductivity is about three order of magnitude larger in the ordered case than in the disorder one. This shows that disorder may have a very strong localization effect.

(0,0) (30,205)[COE]{} (262,205)[CUE]{}

Finally, let us discuss the effects of decoherence. The idea is to view the self-consistent reservoirs as an effective environment acting on a chain of coherent quantum dots. For simplicity, let us assume that all the self-consistent reservoirs are equally coupled to the dots, $1 - |S^{(k)}_{22}|^2 = \lambda \in [0,1]$ for all $k$. (In the simulations, the values of the coupling parameter $\lambda$ are taken in small intervals of width $0.01$.) Then, $\lambda = 0$ corresponds to the completely decoupled situation and $\lambda = 1$ corresponds to the maximally coupled one.

For any fixed value of $\lambda$, we found that the profiles and conductivities are qualitatively similar to those obtained previously. Interestingly, when $\lambda$ increases, we observed that the conductivity increases in the disordered cases while it decreases in the ordered ones. As an illustration, we show in  \[Conductivity Ordered model Coupling\] the results for three different values of $\lambda$ in the quantum disordered chaotic chain. In this case, the profiles (not shown) are similar to those obtained in  \[Fig B profile\] and tend to get closer and closer to linear as the coupling parameter $\lambda$ increases.

(0,0) (125,235)[COE]{}

Concluding Remarks
==================

We have presented a model for charge and heat transport using the Landauer-Büttiker scattering approach. In the linear response regime, we have seen that the Onsager relations hold. Then, assuming that the scattering matrix $S$ does not depend on the energy, we have shown that the transport matrix $L$ is real positive semi-definite and we have characterized the equilibrium and non-equilibrium states of the system in terms of the currents and in terms of the entropy production. In particular, we have shown that our multi-terminal system satisfies the first and second laws of thermodynamics.

In order to obtain an effective quantum version of the Eckmann-Young model with rotating discs, we have imposed the self-consistency condition: $I_i = J_i = 0$ for $i=1,\dots,N$. This condition led to expressions for the temperature and chemical potential profiles, which turn out to be independent of the nature of the particles, as well as to expressions for the electric and heat currents going through the system. Finally, we have presented some numerical results, using random matrix theory, supporting the validity of Ohm and Fourier laws in our model.

Let us point out that we could handle the three physical cases, Maxwell-Boltzmann, Fermi-Dirac and Bose-Einstein, on the same footing because they share the following two features: (i) the specific form of $f$ and (ii) the ratio $\mathcal{R}$, defined in , satisfies $0 < \mathcal{R} <
1$. Basically, we expect that any distribution function $f$ satisfying (i)-(ii) will lead to the same results.

For further investigations, it would be interesting to study the general energy-dependent S-matrix situation and to compare the corresponding predictions (under RMT) with those of the EY-model. One may also analyse the nonlinear transport properties through the full counting statistics (FCS) [@Levitov-Lesovik; @levitov2; @Bagrets; @Pilgram-Jordan; @Pilgram2; @Kindermann; @Pilgram; @Saito1]. More mathematically oriented investigations would be to obtain some rigorous results concerning the global composite scattering matrix $S$ and consequently to prove the existence of the finite limit, $\lim_{N \rightarrow
  \infty}\kappa(N)$, which ensures the validity of Ohm and Fourier laws in our model.

Acknowledgements {#acknowledgements .unnumbered}
================

The author has benefited from numerous helpful discussions with M. Büttiker and his group members, H. Foerster, A. M. Lunde, S. Nigg, C. Petitjean, M. L. Polianski and J. Splettstösser, and also with J.-P. Eckmann, M. Hairer, J. Jacquet, Ph. Jacquod, C. Mejía-Monasterio, M. Moskalets, C.-A. Pillet, L. Rey-Bellet, E. Sukhorukov, P. Wittwer and C. Zbinden, to all of whom he wishes to express his sincere gratitude. This work was partially supported by the Swiss National Science Foundation and the Erwin Schrödinger International Institute for Mathematical Physics.

Appendix A: The Coefficients $\mathbf{C(n)}$ {#The coefficients .unnumbered}
============================================

Here we derive some properties of the coefficients $C(0)$, $C(1)$ and $C(2)$. For this, let $n = 0$, $1$ or $2$. Then\
\
(i) If $x_0 \in (-1,\infty)$, then $C^{\rm MB}(n)
> 0$. Indeed, this follows from $$\label{expl MB}
 C^{\rm MB}(0) = e^{-x_0}~, \hspace{2mm} C^{\rm MB}(1) = (1+x_0)
e^{-x_0}~, \hspace{2mm}C^{\rm MB}(2) = (2+2 x_0+x_0^2~)
e^{-x_0}~.$$ (ii) If $x_0 \in (-\infty,\infty)$, one has $C^{\rm FD}(n) > 0$, and if $x_0 \in (0,\infty)$, one has $C^{\rm BE}(n) > 0$. Indeed, if $x_0 >
0$, then clearly $C^{\rm FD}(n) > 0$ and $C^{\rm BE}(n) > 0$. For the Fermi-Dirac case, assume that $x_0 \leq 0$, then $$\begin{aligned}
C^{\rm FD}(0) &=& \int_{x_0}^{\infty}
  \frac{ e^x}{(e^x + 1)^2} \ dx > 0~,\\
C^{\rm FD}(2) &=& \int_{x_0}^{\infty}
  \frac{x^2 e^x}{(e^x + 1)^2} \ dx > 0~.\end{aligned}$$ Now, $$C^{\rm FD}(1) = \int_{x_0}^{\infty}
  \frac{x e^x}{(e^x + 1)^2} \ dx~.$$ Observe that $C^{\rm FD}(1) > 0$ if $x_0 = 0$ and $C^{\rm FD}(1) \rightarrow 0$ as $x_0 \rightarrow -\infty$. Therefore, the property $C^{\rm FD}(1) > 0$, for all $x_0 \leq 0$, follows from $$\frac{d}{dx_0} C^{\rm FD}(1) = -\frac{x_0 e^{x_0}}{(e^{x_0} + 1)^2} \geq 0~.$$

Appendix B: The Ratio $\mathbf{\mathcal{R}}$ {#appendix-b-the-ratio-mathbfmathcalr .unnumbered}
============================================

Here, we show that for all $x_0$: $$\label{Lemma L0}
\mathcal{R} \equiv \frac{Q_1}{\sqrt{Q_2}} =
\frac{C(1)}{\sqrt{C(0) \cdot C(2)}} \in (0,1)~.$$ Since $Q_1 > 0$ and $Q_2 > 0$, one has $\mathcal{R} > 0$, so it only remains to show that $\mathcal{R} < 1$.\
\
(i) Let us first consider the Maxwell-Boltzmann case. Recalling the explicit expressions , one obtains ($x_0 \in (-1,\infty)$) $$\mathcal{R}^2 = \frac{1 + 2 x_0 + x_0^2}{2 + 2 x_0 + x_0^2} < 1~.$$ To show that $\mathcal{R} < 1$ in the other two cases, it is convenient to introduce the real functional space $L^2([a,\infty),dx)$ with scalar product $(\cdot, \cdot)_{a}$ and norm $\| \cdot \|_{a}$: $$(f_1, f_2)_{a} = \int_{a}^{\infty} f_1(x) \ f_2(x) \ dx~,
  \hspace{10mm} \|
  f\|^2_{a} = (f, f)_{a}~.$$ (ii) For the Bose-Einstein case, one has ($x_0 > 0$) $$\mathcal{R} = \frac{\| \sqrt{x} \ g \|_{x_0}^2}{\|g \|_{x_0} \cdot
    \|x \ g \|_{x_0}}~, \hspace{5mm} \mbox{with} \hspace{5mm} g(x) = \frac{e^{x/2}}{|e^x - 1|}~.$$ Hence, it is sufficient to show that $$\| \sqrt{x} \ g \|_{x_0}^2 < \|g \|_{x_0} \cdot \| x \ g
\|_{x_0}~.$$ We have $$\| \sqrt{x} \ g \|_{x_0}^2 = (\sqrt{x} \ g, \sqrt{x} \ g)_{x_0} = (g, x \ g)_{x_0}
< \| g \|_{x_0} \cdot \| x \ g \|_{x_0}~,$$ where we have used the Cauchy-Schwarz inequality.

\(iii) Finally, let us consider the Fermi-Dirac case. If $x_0 > 0$, then we proceed as in (ii). If $x_0 \leq  0$, then we observe that $$\mathcal{R} \leq \frac{\| \sqrt{x} \ g \|_{0}^2}{\|g \|_{0} \cdot
    \|x \ g \|_{0}}~, \hspace{5mm} \mbox{with} \hspace{5mm} g(x) = \frac{e^{x/2}}{|e^x + 1|}~,$$ and again proceed as in point (ii).

Appendix C: The Entropy Production {#Appendix Reminder .unnumbered}
==================================

In this appendix, we derive the expression for the entropy production rate $\sigma_s$. However, let us first recall how a real positive semi-definite matrix is defined. Let $\{e_j\} \subset \mathbb{R}^{2N+4}$ denote the canonical basis in which $(e_i, L e_j) = L_{ij}$, where $(\cdot, \cdot)$ denotes the usual scalar product in $\mathbb{R}^{2N+4}$. Then the matrix $L$ is said to be real positive semi-definite if for any vector $V
\in \mathbb{R}^{2N+4}$ one has $$\sigma_s = (V, L V)  = \sum_{i} \sum_{j} L_{ij} V_i V_j \geq 0~.$$ Let us decompose the matrix elements $L_{ij}$ into a symmetric and antisymmetric part: $$L_{ij} = L^s_{ij} + L^a_{ij}~,$$ where $$L^s_{ij} = \frac{1}{2} (L_{ij} + L_{ji}) = L^s_{ji} \hspace{5mm} \mbox{and}
\hspace{5mm} L^a_{ij} = \frac{1}{2} (L_{ij} - L_{ji}) = -L^a_{ji}~.$$ Then $$\label{Entropy}
\sigma_s = \sum_{i} \sum_{j} L^s_{ij} V_i V_j~.$$ This shows that only the symmetric part of the transport matrix $L$ may contribute to the entropy production. Observe that the symmetric transport coefficients, $L^s_{ij}$, also satisfy the relations and . Since we shall only make use of these properties, we see that we can assume that the matrix $L$ is symmetric.

The first idea is to decompose the sum into four sums, one associated to each block appearing in the matrix $L$. To simplify the notation, let us consider that the chain contains $N-2$ quantum dots and then number as $1$ and $N$ the left and right reservoirs, respectively. We shall then write $X_1, \dots, X_N$ for the first $N$ components of $V$ and $Y_1, \dots, Y_N$ for its last $N$ components. One may interpret this decomposition as follows: $X_i = \delta \mu_i/e$ and $Y_i = \delta
T_i / T$. Recalling that $L_{ij} = L_{ji}$, one has $$\sigma_s  = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i X_j +
\sum_{i=1}^{N} \sum_{j=1}^{N} (L^{(1)}_{ij} + L^{(1)}_{ji}) X_i Y_j +
\sum_{i=1}^{N} \sum_{j=1}^{N} L^{(2)}_{ij} Y_i Y_j~.$$ Now, using the relations $L^{(1)}_{ij} = Q_1L^{(0)}_{ij}  = L_{ji}^{(1)}$ and $L^{(2)}_{ij} = Q_2
L^{(0)}_{ij}$, one obtains $$\sigma_s = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} \left[X_i X_j +
2 Q_1 X_i Y_j + Q_2 Y_i Y_j\right]~.$$ Let $Z_j = \sqrt{Q_2} Y_j$ and $\mathcal{R} =
Q_1/\sqrt{Q_2}$. Then $$\sigma_s = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} \left[X_i X_j +
2 \mathcal{R} X_i Z_j + Z_i Z_j\right]~.$$ Let us first consider the term: $$\mathcal{T}_1 = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i X_j~.$$ We write $$\mathcal{T}_1 = \sum_{i=1}^{N} L^{(0)}_{ii} X_i^2 + \sum_{i=1}^{N} \sum_{j=1}^{N}
L^{(0)}_{ij} X_i X_j (1 - \delta_{ij})~.$$ Since $\sum_{j} L^{(0)}_{ij} = 0$, one has $L^{(0)}_{ii} = - \sum_{j}
L^{(0)}_{ij} (1-\delta_{ij})$ and therefore $$\mathcal{T}_1 = -\sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i^2 (1-\delta_{ij})+ \sum_{i=1}^{N} \sum_{j=1}^{N}
L^{(0)}_{ij} X_i X_j (1 - \delta_{ij})~.$$ Now comes a trick: $$\sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i^2 (1-\delta_{ij}) = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_j^2 (1-\delta_{ij})~.$$ Hence $$\mathcal{T}_1 = -\sum_{\scriptsize\begin{array}{c}
    i,j = 1\\ i < j \end{array}}^{N} L^{(0)}_{ij} \underbrace{(X_i^2
    -2 X_i X_j + X_j^2)}_{(X_i- X_j)^2}~.\label{Equ XX}$$ This relation is well known (see [@Buttiker-Symmetry]). Next we consider the second term: $$\mathcal{T}_2 = 2 \mathcal{R} \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i Z_j
= 2 \mathcal{R} \sum_{i=1}^{N} L^{(0)}_{ii} X_i Z_i + 2 \mathcal{R} \sum_{i=1}^{N} \sum_{j=1}^{N}
L^{(0)}_{ij} X_i Z_j (1 - \delta_{ij})~.$$ Using the relation $L^{(0)}_{ii} = - \sum_{j}
L^{(0)}_{ij} (1-\delta_{ij})$ and the identities $$\sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i Z_i (1 -
\delta_{ij}) = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_j Z_j (1 - \delta_{ij})~,$$ and $$\sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_i Z_j (1 -
\delta_{ij}) = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} X_j Z_i (1 - \delta_{ij})~,$$ one can write $$\mathcal{T}_2 = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} \mathcal{R} \left(X_i
  Z_j + X_j Z_i - X_i Z_i -  X_j Z_j\right) (1 - \delta_{ij})~.$$ Writing the term in $Z_i Z_j$ in a similar manner as $\mathcal{T}_1$, one finally obtains $$\label{expr sigmas}
\sigma_s = \sum_{\scriptsize\begin{array}{c}
    i,j = 1\\ i < j \end{array}}^{N} (-L^{(0)}_{ij}) \ I_{ij}~,$$ where $$I_{ij} = (X_i - X_j)^2
  + (Z_i - Z_j)^2 -  2 \mathcal{R} C_{ij}~, \hspace{5mm} C_{ij} = X_i Z_j + X_j Z_i - X_i Z_i - X_j Z_j~.$$ This is the relation used in  \[The entropy production\] to prove that $\sigma_s \geq 0$.

Appendix D: Equilibrium and Non-Equilibrium States {#appendix-d-equilibrium-and-non-equilibrium-states .unnumbered}
==================================================

In this appendix, we prove the equivalent relations presented in  \[Equilibrium and non-equilibrium states\].

To obtain the first one, $$\{\mbox{System is at equilibrium}\} \iff \{I_i = 0 \hspace{3mm}
\mbox{and} \hspace{3mm} J_i = 0, \hspace{3mm} \forall i\}~,$$ let us write the no-current condition ($i\in\{\L, \R, 1, \dots, N\}$): $$\begin{aligned}
I_i &=& \sum_j L_{ij}^{(0)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(1)}
\, \frac{\delta T_j}{T} = 0~,\label{cond 1 I}\\
J_i &=& \sum_j L_{ij}^{(1)} \, \frac{\delta \mu_j}{e} + L_{ij}^{(2)}
\, \frac{\delta T_j}{T} = 0~.\label{cond 1 J}\end{aligned}$$ It is convenient to define $X, Y \in \mathbb{R}^{N+2}$ by $$X_j = \frac{\delta \mu_j}{e} \hspace{5mm} \mbox{and} \hspace{5mm} Y_j
= \frac{\delta T_j}{T}~.$$ Then, using the relations $L_{ij}^{(1)} = Q_1 L_{ij}^{(0)}$ and $L_{ij}^{(2)} = Q_2
L_{ij}^{(0)}$, and setting $Q = Q_2 / Q_1 $, one sees that the conditions – can be rewritten as follows: $$\begin{aligned}
\label{XY}
L^{(0)} (X + Q_1  Y) &=& 0~,\label{XY1} \\
L^{(0)} (X + Q  Y)&=& 0~.\label{XY2}\end{aligned}$$ Since $\mathcal{R} \not = 1$, it follows that $Q_1
\not = Q$. Assume now that all the currents vanish. Then the equations – imply that $$L^{(0)} X = 0 \hspace{5mm} \mbox{and} \hspace{5mm} L^{(0)} Y = 0~.$$ Since $\sum_j L^{(0)}_{ij} = 0$, the equilibrium situation, $X_i = \delta \mu/e$ and $Y_i
= \delta T/T$, is a solution. Observe now that $(-L^{(0)})$ may be interpreted as the generator of an irreducible Markov chain and, consequently, the equation $L^{(0)} V = 0$ has a unique (normalised) solution. Hence the system is at equilibrium. Finally, the reverse implication is immediate since the equilibrium situation satisfies the equations –.

We next turn to the equivalences in terms of the entropy production: $$\left\{\hspace{-1mm}\begin{array}{cc}\mbox{System is at}\\\mbox{equilibrium}\end{array}\hspace{-1mm}\right\} \iff \sigma_s = 0 \hspace{5mm} \mbox{and}
\hspace{5mm} \left\{\hspace{-1mm}\begin{array}{cc}\mbox{System is}\\\mbox{out of
  equilibrium} \end{array}\hspace{-1mm}\right\}\iff \sigma_s > 0~.$$ Observe that, since $\sigma_s
 \geq 0$, the two above statements are contrapositive equivalences. Therefore, it is sufficient to prove the first. For this we use the expression for $\sigma_s$.

\[$\Rightarrow$\] Assume the system is at equilibrium. Then, $X_1 = \dots = X_{N}$ and $Z_{1} = \dots =
 Z_{N}$, and therefore $I_{ij} = 0$ for all $i,j$, from which it follows that $\sigma_s = 0$.

\[$\Leftarrow$\] Conversely, assume that $\sigma_s = 0$. Then, since $(-L^{(0)}_{ij}) < 0$ and $I_{ij} \geq 0$ for all $i < j$, one has $I_{ij} = 0$ for all $i < j$. Therefore, by (\[In1a\]), one deduces that $X_1 = \dots = X_{N}$ and $Z_{1} = \dots =
 Z_{N}$, the system is at equilibrium.

Appendix E: The Matrix $\mathbf{L^{(0)}_{\rm C}}$ {#appendix-e-the-matrix-mathbfl0_rm-c .unnumbered}
=================================================

Here we prove that the matrix $L^{(0)}_{\rm C}$ is real positive definite and consequently that $L^{(0)}_{\rm C}$ has positive determinant and is therefore invertible.

Let $W \in \mathbb{R}^{N}$, $W \not = 0$, and $$\mathcal{T} \equiv (W, L^{(0)}_{\rm C} W)  = \sum_{i=1}^{N} \sum_{j=1}^{N} L^{(0)}_{ij} W_i W_j~.$$ The idea is to proceed as in the proof of $\sigma_s \geq 0$, but with the identity $$L^{(0)}_{ii} = - L^{(0)}_{i\L} - L^{(0)}_{i\R} - \sum_{j=1}^{N}
L^{(0)}_{ij} (1-\delta_{ij})~.$$ We obtain $$\mathcal{T} = -\sum_{\scriptsize\begin{array}{c}
    i,j = 1\\ i < j \end{array}}^{N} L^{(0)}_{ij} (W_i- W_j)^2 -
    \sum_{i=1}^{N} (L^{(0)}_{i\L} + L^{(0)}_{i\R}) W_i^2 > 0~.$$ This shows that $L^{(0)}_{\rm C}$ is real positive definite.

Now, let $\lambda_1,
\dots, \lambda_N$ denote the eigenvalues of $L^{(0)}_{\rm C}$, which might be complex since $L^{(0)}_{\rm C}$ is not symmetric in general.

-   Assume that $\lambda$ is a real eigenvalue of $L^{(0)}_{\rm C}$ and let $v \in
      \mathbb{R}^{N}$ be a normalised associated eigenvector. Then $$\lambda = \lambda \| v\| = (v, \lambda v) = (v,  L^{(0)}_{\rm C}
    v) > 0~,$$ where we have used that $L^{(0)}_{\rm C}$ is real positive definite.

-   Assume next that $\lambda$ is a complex eigenvalue: $\lambda = a + i b$, with $a,b \in \mathbb{R}$ and $b \not = 0$. Then, since $L^{(0)}_{\rm C}$ is a real matrix, both $\lambda$ and its complex conjugate $\overline{\lambda}$ are eigenvalues of $L^{(0)}_{\rm C}$.

In summary, we have shown that the eigenvalues of $L^{(0)}_{\rm C}$ are either positives or come in pairs $(\lambda, \overline{\lambda})$ with $\lambda \not = 0$. Therefore $$\det(L^{(0)}_{\rm C}) = \prod_i \lambda_i > 0~.$$

Appendix F: The Profiles {#appendix-f-the-profiles .unnumbered}
========================

Here we explain the details in the derivation of the profiles –. From , one sees that the profiles are obtained by solving the generic equation: $$L^{(0)}_{\rm C} X = - \sum_{\ell = \L, \R}  D_\ell X_\ell~.$$ From Appendix E, the matrix $L^{(0)}_{\rm C}$ can be inverted to give $$X = - \sum_{\ell = \L, \R} [L^{(0)}_{\rm C}]^{-1} D_\ell X_\ell~.$$ Recalling that $(D_{\ell})_j = L^{(0)}_{j \ell}$, one has in components: $$X_i = - \sum_{j = 1}^{N} ([L^{(0)}_{\rm C}]^{-1})_{ij} \left[L^{(0)}_{j \L}
  X_\L + L^{(0)}_{j \R}
  X_\R \right]~.$$ Now comes a trick: $$L^{(0)}_{j \L} = - L^{(0)}_{j \R} - \sum_{k = 1}^{N} L^{(0)}_{j k}~.$$ With this identity, one has $$\begin{aligned}
X_i &=& \sum_{k = 1}^{N} \sum_{j = 1}^{N} ([L^{(0)}_{\rm C}]^{-1})_{ij}
L^{(0)}_{j k} X_\L - \sum_{j = 1}^{N} ([L^{(0)}_{\rm C}]^{-1})_{ij}
L^{(0)}_{j \R} (X_\R - X_\L)\nonumber\\
&=& X_\L - \sum_{j = 1}^{N} ([L^{(0)}_{\rm C}]^{-1})_{ij}
L^{(0)}_{j \R} (X_\R - X_\L)~.\end{aligned}$$ Observe that $(L^{(0)}_{C})_{ij} =
  L^{(0)}_{ij} = e^2 C(0)/h \ \Gamma_{ij}$ and $L^{(0)}_{j\R} = - e^2 C(0)/h \ t_{j\R}$. Hence, the constant $e^2 C(0)/h$ involved in $([L^{(0)}_{\rm C}]^{-1})_{ij}$ and in $L^{(0)}_{j\R}$ simplifies, and we are left with profiles independent of the distribution $f$ describing the reservoirs. The temperature and chemical potential profiles are obtained by setting $X_i = \delta T_i / T$, with $T_i =
T + \delta T_i$, and $X_i = \delta \mu_i / e$, with $\mu_i =
\mu + \delta \mu_i$, respectively.

Appendix G: An Expression for $\mathbf{A_i}$ {#Appendix Reminder .unnumbered}
============================================

In this appendix, we derive the alternative expression for the coefficient $A_i$. From , one has $$A_i = \sum_{j = 1}^{N} (\Gamma_{\rm C}^{-1})_{ij} \ t_{j \R}~.$$ Now, we write $$\label{equ gammacinv}
(\Gamma_{\rm C}^{-1})_{ij} = \frac{1}{\det\left(\Gamma_{\rm C}\right)} (-1)^{i+j} \det\left(\Gamma_{\rm C}(j,i)\right)~,$$ so that $$A_i = \frac{1}{\det\left(\Gamma_{\rm C}\right)} \sum_{j = 1}^{N} (-1)^{i+j}
\det\left(\Gamma_{\rm C}(j,i)\right) t_{j \R}~.$$ It remains to expand $\det\left(\Gamma_{\rm C}\right)$. One has $$\det\left(\Gamma_{\rm C}\right) =  \sum_{j = 1}^{N} (-1)^{i+j} \det\left(\Gamma_{\rm C}(j,i)\right) \Gamma_{ji}~.$$ Using the relation $$\Gamma_{ji} = - \Gamma_{j \L} - \Gamma_{j \R} - \sum_{\scriptsize\begin{array}{c}
    k = 1\\ k \not
= i\end{array}}^{N} \Gamma_{j k}~,$$ one obtains $$\det\left(\Gamma_{\rm C}\right) = \sum_{j = 1}^{N} (-1)^{i+j}
\det\left(\Gamma_{\rm C}(j,i)\right) [t_{j \L} + t_{j \R}] - R_i~,$$ where $$R_i = \sum_{j = 1}^{N} \sum_{\scriptsize\begin{array}{c}
    k = 1\\ k \not
= i\end{array}}^{N} (-1)^{i+j}
\det\left(\Gamma_{\rm C}(j,i)\right) \Gamma_{j k}~.$$ Finally, using , one obtains $$R_i = \det\left(\Gamma_{\rm C}\right) \sum_{\scriptsize\begin{array}{c}
    k = 1\\ k \not
= i\end{array}}^{N} \sum_{j = 1}^{N}  \left(\Gamma_{\rm
    C}^{-1}\right)_{ij}  \Gamma_{j k} = \det\left(\Gamma_{\rm
    C}\right) \sum_{\scriptsize\begin{array}{c}
    k = 1\\ k \not
= i\end{array}}^{N} \delta_{ik} = 0~.$$

Appendix H: The Global Scattering Matrix $\mathbf{S_N}$ {#General case S .unnumbered}
=======================================================

In this appendix, we explain how to build the scattering matrix $S$ of the system made of $N$ quantum dots from the $N$ individual scattering matrices $S^{(1)},\dots,S^{(N)}$, where $S^{(k)}$ is the $(M + M_k + M) \times (M
+ M_k + M)$ unitary scattering matrix associated with the $k$-th quantum dot. Here, we have set $M_\L = M_\R = M$.

Let us start with the case $N=2$. Let $S^{(1)}$ and $S^{(2)}$ denote the scattering matrices of the first and second quantum dot, respectively, and let $S_2$ denote the corresponding scattering matrix for the system made of two quantum dots (see  \[Figure S\]).

From  \[Figure S\], one easily sees that these matrices are related as follows: $$\label{N equal to 1}
S^{(1)} {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}x_\L \\ x_1 \\ z \\ \end{array}\hspace{-2mm}\right)} =
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}y_\L \\ y_1 \\ z' \\ \end{array}\hspace{-2mm}\right)}~, \hspace{5mm} S^{(2)} {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}z' \\ x_2 \\ x_\R \\ \end{array}\hspace{-2mm}\right)} =
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}z \\ y_2 \\ y_\R \\ \end{array}\hspace{-2mm}\right)}~, \hspace{5mm} S_{2}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}x_\L \\ x_1 \\ x_2 \\ x_\R \\ \end{array}\hspace{-2mm}\right)} =
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}y_\L \\ y_1 \\ y_2 \\ y_\R \\ \end{array}\hspace{-2mm}\right)}~.$$ Introducing $$r^{(1)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(1)}_{11} & S^{(1)}_{12} \\ S^{(1)}_{21}
  & S^{(1)}_{22}\\ \end{array}\hspace{-2mm}\right)}~, r'^{(1)} = S^{(1)}_{33}~, t^{(1)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(1)}_{31} & S^{(1)}_{32} \\ \end{array}\hspace{-2mm}\right)}~,  t'^{(1)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(1)}_{13} \\ S^{(1)}_{23} \\ \end{array}\hspace{-2mm}\right)}$$ and $$x_{\L 1} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} x_\L \\ x_1 \\ \end{array}\hspace{-2mm}\right)}~,
y_{\L 1} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} y_\L \\ y_1 \\ \end{array}\hspace{-2mm}\right)}$$ one can rewrite the first relation in as follows: $$\label{Equ S1}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} r^{(1)} & t'^{(1)} \\ t^{(1)} & r'^{(1)}\\ \end{array}\hspace{-2mm}\right)}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} x_{\L 1} \\ z \\ \end{array}\hspace{-2mm}\right)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} y_{\L 1} \\ z' \\ \end{array}\hspace{-2mm}\right)}~.$$ Similarly, setting $$r^{(2)} = S^{(2)}_{11}~,   r'^{(2)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(2)}_{22} & S^{(2)}_{23} \\ S^{(2)}_{32}
  & S^{(2)}_{33}\\ \end{array}\hspace{-2mm}\right)}~, t^{(2)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(2)}_{21} \\ S^{(2)}_{31} \\ \end{array}\hspace{-2mm}\right)}~,  t'^{(2)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(2)}_{12} & S^{(2)}_{13} \\ \end{array}\hspace{-2mm}\right)}$$ and $$x_{2 \R} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} x_2 \\ x_\R \\ \end{array}\hspace{-2mm}\right)}~,
y_{2 \R} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} y_2 \\ y_\R \\ \end{array}\hspace{-2mm}\right)}$$ one has for the second relation in : $$\label{Equ S2}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} r^{(2)} & t'^{(2)} \\ t^{(2)} & r'^{(2)}\\ \end{array}\hspace{-2mm}\right)}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} z' \\  x_{2 \R}\\ \end{array}\hspace{-2mm}\right)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} z \\  y_{2 \R} \\ \end{array}\hspace{-2mm}\right)}~.$$ Finally, we write the third relation in as follows: $$\label{Equ S12}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} r_{2} & t'_{2} \\ t_{2} & r'_{2}\\ \end{array}\hspace{-2mm}\right)}
{\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} x_{\L 1} \\  x_{2 \R}\\ \end{array}\hspace{-2mm}\right)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}
  y_{\L 1} \\  y_{2 \R} \\ \end{array}\hspace{-2mm}\right)}~.$$ Solving the equations  and for $y_{\L 1}$ and $y_{2 \R}$, and recalling , one obtains $$\begin{aligned}
r_{2} &=& r^{(1)} + t'^{(1)} [1 - r^{(2)} r'^{(1)}]^{-1} r^{(2)}
t^{(1)}~,\\
r'_{2} &=& r'^{(2)} + t^{(2)} [1 - r'^{(1)} r^{(2)}]^{-1} r'^{(1)}
t'^{(2)}~,\\
t_{2} &=& t^{(2)} [1 - r'^{(1)} r^{(2)}]^{-1} t^{(1)}~,\\
t'_{2} &=& t'^{(1)} [1 - r^{(2)} r'^{(1)}]^{-1} t'^{(2)}~.\end{aligned}$$ Let us now generalise these relations to the $N > 2$ quantum dots situation. For this we denote by $S^{(1)},\dots,S^{(N)}$ the $N$ individual scattering matrices and by $S_{N}$ the composite S-matrix associated to $N$ quantum dots. We write $$S_{N} {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}x_\L \\ x_1 \\ : \\ x_N \\x_\R
\\ \end{array}\hspace{-2mm}\right)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc}y_\L \\ y_1 \\ : \\ y_N \\ y_\R
\\ \end{array}\hspace{-2mm}\right)}~, \hspace{10mm} S_{N} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} r_{N} & t'_{N} \\ t_{N} & r'_{N}\\ \end{array}\hspace{-2mm}\right)}~.$$ Then $$\begin{aligned}
r_{N} &=& r^{(12\dots N-1)} + t'^{(12\dots N-1)} [1 - r^{(N)} r'^{(12\dots N-1)}]^{-1} r^{(N)}
t^{(12\dots N-1)}~,\label{Eq rN}\\
r'_{N} &=& r'^{(N)} + t^{(N)} [1 - r'^{(12\dots N-1)} r^{(N)}]^{-1} r'^{(12\dots N-1)}
t'^{(N)}~,\\
t_{N} &=& t^{(N)} [1 - r'^{(12\dots N-1)} r^{(N)}]^{-1} t^{(12\dots N-1)}~,\\
t'_{N} &=& t'^{(12\dots N-1)} [1 - r^{(N)} r'^{(12\dots N-1)}]^{-1} t'^{(N)}~,\label{Eq tpN}\end{aligned}$$ where $$r^{(N)} = S^{(N)}_{11}~, \; r'^{(N)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(N)}_{22} & S^{(N)}_{23} \\ S^{(N)}_{32}
  & S^{(N)}_{33}\\ \end{array}\hspace{-2mm}\right)}~, \; t^{(N)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(N)}_{21} \\ S^{(N)}_{31} \\ \end{array}\hspace{-2mm}\right)}~,  \;t'^{(N)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} S^{(N)}_{12} & S^{(N)}_{13} \\ \end{array}\hspace{-2mm}\right)}~.$$ and $$r^{(12\dots N-1)} = {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} (S_{N-1})_{11} & \dots &
  (S_{N-1})_{1N} \\ : & & : \\
  (S_{N-1})_{N1} & \dots & (S_{N-1})_{NN}\\ \end{array}\hspace{-2mm}\right)}~, \; r'^{(12\dots N-1)} =
  (S_{N-1})_{N+1 \ N+1}~,$$ $$t^{(12\dots N-1)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} (S_{N-1})_{N+1 \ 1} & \dots &  (S_{N-1})_{N+1 \ N}\\ \end{array}\hspace{-2mm}\right)}~, \; t'^{(12\dots N-1)} =
  {\left(\hspace{-2mm}\begin{array}{cccccccccccccccccc} (S_{N-1})_{1 \ N+1} \\ : \\ (S_{N-1})_{N \ N+1} \\ \end{array}\hspace{-2mm}\right)}~.$$ Therefore, by working recursively, one can build the global scattering matrix $S = S_{N}$ from the local scattering matrices $S^{(1)},\dots,S^{(N)}$. Note that, by construction, if $S^{(1)},\dots,S^{(N)}$ are unitary, then $S$ will also be unitary.\
\
*Remark.* If one considers that the terminals have only one channel: $M_i = 1$ for all $i$, then $r^{(N)}$ and $r'^{(12\dots N)}$ become complex numbers and we can introduce the following complex number: $$c_N = [1 - r^{(N)} r'^{(12\dots N-1)}]^{-1}~.$$ Then one can rewrite the equations – as follows: $$\begin{aligned}
r_{N} &=& r^{(12\dots N-1)} + c_N \ r^{(N)} \ t'^{(12\dots N-1)} t^{(12\dots N-1)}~,\label{Eq rN2}\\
r'_{N} &=& r'^{(N)} + c_N \ r'^{(12\dots N-1)} \ t^{(N)} t'^{(N)}~,\\
t_{N} &=& c_N \ t^{(N)} t^{(12\dots N-1)}~,\\
t'_{N} &=& c_N \ t'^{(12\dots N-1)} t'^{(N)}~.\label{Eq tpN2}\end{aligned}$$ These are the equations that are used in the numerical simulations in  \[Numerical analysis\].

[^1]: This way of modelling the discs in the scattering approach was suggested to me by M. Büttiker.
