M Lavagna

V Talbo

T Q Duong

Adeline Crépieux

Level anticrossing effect in single-level or multilevel double quantum dots: Electrical conductance, zero-frequency charge susceptibility and Seebeck coefficient

HAL is

I. INTRODUCTION

The study of double quantum dot (DQD) has been the focus of an increasing number of research works in the last years both theoretically and experimentally. One of the main reasons explaining these developments are that DQDs are promising candidates to build quantum bits of spin, i.e. spin qubits 1,2 , which surpass charge qubits because of their longer coherence time 3 . Moreover, a DQD is a readily tunable system [4][5][6] : the number of electrons in each of the two dots can be controlled by varying gate voltages located at proximity. Proposals to probe and drive the spin and charge states in DQDs have also been made [7][8][9] . Note that the second dot is used under certain circumstances to control the spin state in the first dot benefiting from the Pauli spin blockade effect 10 . Initially, experimental DQDs were built from GaAs heterostructures [11][12][13][14][15][16][17] but one has witnessed in the last five years the development of studies in Si-based DQDs which have the advantage over the former one to present a longer spin coherence time [18][19][20][21][22][23][24][25][26][27][28][29] . The serialcoupling of the two dots is by far the geometry that has been the most studied. However, the case of two parallel coupled dots is interesting too since it may give rise to interference effects [30][31][32] or other specific effects [33][34][35][36][37][38][39][40] . In this paper, both geometries of serial-coupled and parallelcoupled DQDs are considered. From the theoretical side, the electrical transport properties in DQDs have been widely studied by using various methods going from scattering matrix theory [41][42][43] for noninteracting case to Master or Bloch equa-tion approaches [44][45][46][47][48][49][50] and nonequilibrium Green function methods [51][52][53][54][55][56] for more general cases. A classical theory has also been developed along which the DQD is modeled as a network of resistors and capacitors which mimic the tunnel and electrostatic couplings between dots and leads [4][5][6] . The obtained results for the overall evolution of the linear electrical conductance as a function of gate voltages are as follows: (i) at weak interdot coupling, conductance peaks are observed at the nodes of a square lattice in the (ε 1 , ε 2 ) phase space, where ε 1 and ε 2 are the level energies of each of the two dots, delimiting domains with different integer occupancies in the charge stability diagram, (ii) at intermediate interdot coupling, the nodes separate into pairs of triple points, corresponding to the deformation of the square lattice into a honeycomb lattice in the (ε 1 , ε 2 ) phase space, and (iii) at strong interdot coupling, the triple point separation reaches its maximum and the DQD behaves as a single dot with an occupancy N 1 + N 2 , where N 1 and N 2 are the average rate of electronic occupancy in the dots 1 and 2 respectively.

Charge susceptibility and thermoelectrical properties in DQDs have aroused much less attention than electrical conductance while the growing-up activity in both gate-reflectometry experiments 23,25,57 and thermopower measurements [58][59][60][61][62] in spin qubits provides a strong motivation to intensify the theoretical efforts in that direction. From the theoretical side, one essentially refers to the works of Refs. 63-65 where the charge susceptibility of a DQD is discussed in the context of mesoscopic admittance in either the noninteracting case or the inter-acting case where it is affected by the Pauli spin blockade effect. These works are in line with the seminal works of Büttiker and collaborators 66,67 which demonstrate the importance of taking the charge susceptibility contribution into account in the mesoscopic capacitance of the system, in addition to the standard geometrical capacitances of the dots. The charge susceptibility reflects the ability of the system when it is open, to adjust the average rate of electronic occupancy in the dots to the external excitation brought by the gate voltage V g . It is the relevant physical quantity behind reflectometry experiments carried out in spin qubit systems and plays a key role in the context of manipulation, coupling and readout of qubits. Theoretical works devoted to the thermoelectrical properties of DQDs have highlighted specific features such as the increase of the figure of merit when the interdot coupling is reduced 68 , the decrease of the efficiency at maximum power in the presence of Coulomb interactions in the dots 69 . Some other works have focused on the effects of electron-hole symmetry 70,71 and quantum interferences 72,73 on thermopower. Moreover, it has been shown that a DQD constitutes a minimal thermoelectric generator 74 and the issue related to energy harvesting is a central one [75][76][77] . All these results demonstrate the need to develop further theoretical studies on charge susceptibility and thermopower in DQDs. In this paper, we simultaneously discuss electrical transport, charge susceptibility and thermoelectrical properties of a DQD connected to two reservoirs (leads) of electrons, whether the two dots are coupled in series or in parallel, and contain a single energy level or multiple energy levels. The calculations are performed by using the nonequilibrium Green function technique. We study the evolution of the electrical conductance, Seebeck coefficient, and zero-frequency charge susceptibility with the gate voltages applied to the two dots as well as the stability phase diagram, in various regimes going from weak to strong interdot coupling. The results are valid at any temperature, lead-dot and interdot couplings, bias and gate voltages. They are in qualitative agreement with experiments, with notably the succession of sign changes of the Seebeck coefficient when varying the gate voltages applied to the dots. The plan of the paper is the following. In Sec. II A, we present the hamiltonian describing the DQD system. We adopt a unified presentation which enables to describe the case when the dots are coupled in series as well as the case in parallel. In Sec. II B, we give the analytical expressions for the electrical current which allows us to derive the linear electrical conductance and the Seebeck coefficient. In Sec. II C, we derive the expression for the zero-frequency charge susceptibility related to the occupancies of the dots. The numerical results obtained for conductance, Seebeck coefficient, occupancies of the dots and zero-frequency charge susceptibility are presented in Sec. III in the case when the dots are coupled in series, and in Sec. IV when the dots are coupled in parallel. Section V is a conclusion. Details about the determina- tion of the Green functions in the dots are reported in Appendix A, whereas Appendix B presents the diagonalization of the hamiltonian describing the DQD.

II. MODEL AND EXPRESSIONS FOR ELECTRICAL CURRENT AND CHARGE SUSCEPTIBILITY

A. Hamiltonian

We consider two quantum dots 1 and 2 coupled through a tunnel barrier and connected to two metallic leads as depicted in Fig. 1. Each of the two quantum dot contains N ε and M ε discrete levels of energies ε 1n and

ε 2m respec- tively with n ∈ [0, N ε -1] and m ∈ [0, M ε -1].
H dots = i=1,2 n∈i ε in d † in d in + n∈1 m∈2 V 1n,2m d † 2m d 1n + h.c. (1)
H leads = α=L,R k∈α ε αk c † αk c αk (2)
H hop = α=L,R k∈α i=1,2 n∈i V in,αk c † αk d in + h.c. (3)
in which d † in and d in are the creation and annihilation operators of one electron in the dot i with energy ε in , the index i taking the values 1 and 2. c † αk and c αk are the creation and annihilation operators of one electron in the lead α with momentum k and energy ε αk , the index α taking the values L and R. V 1n,2m and V in,αk are the hopping matrix elements between the states |1n and |2m in the dots and those between the states |in in the dot i and |αk in the lead α. We have V * 1n,2m = V 2m,1n and V * in,αk = V αk,in . The abbreviation h.c. stands for hermitian conjugate. We want to emphasize the very general feature of the expression we introduce for the hamiltonian. It allows one to describe all the possible assemblies of two dots in a unified way: the serial assembly corresponding to the case where V 2m,Lk = V 1n,Rk = 0 for any index n, m or k (depicted in Fig. 1(a)), and the parallel one corresponding to the case where V in,αk = 0 for any index i, n, α or k (depicted in Fig. 1(b)). The spin degree of freedom can be added without any difficulty, which is essential if one wants to describe spin qubits or magnetic leads for instance. However, it will not be included in this paper since we restrict our study to a noninteracting DQD system connected to nonmagnetic leads.

B. Electrical current

In this section we give the expression for the electrical current assuming that the noninteracting DQD is in a steady state. For that we start from the current operator from the lead α defined as I α (t) = -ed N α (t)/dt where N α (t) = k∈α c † αk (t) c αk (t) is the average number of electrons in the lead α. In the steady state the derivative with respect to the time variable is given by d N α (t)/dt = [ N α (t), H]/i . In the limit of wide flat band and energy-independent hopping integrals, one obtains

I α = e h ∞ -∞ Tr Γ α G r (ε) Γ α G a (ε) × f α (ε) -f α (ε) dε (4)
where

α = R for α = L and α = L for α = R, f α (ε) is the Fermi-Dirac distribution function in the lead α, G r(a) (ε)
is the retarded (advanced) Green function matrix in the dots defined as

G r(a) (ε) = G r(a) 11 (ε) G r(a) 12 (ε) G r(a) 21 (ε) G r(a)
22 (ε) (5) and Γ α is the dot-lead coupling matrix defined as

Γ α = 2πρ α |V 1α | 2 V * 1α V 2α V 1α V * 2α |V 2α | 2 (6)
with ρ α , the density of states in the lead α. Tr[ ] is the trace of the matrix. Note that the matrix Γ α is diagonal in the case of a DQD coupled in series and nondiagonal in the case in parallel. Equation (4) corresponds to the Landauer formula for the electrical current with a transmission coefficient equal to

T αα (ε) = Tr Γ α G r (ε)Γ α G a (ε) (7)
Appendix A gives the details of these calculations performed in the framework of the nonequilibrium Green function technique [78][79][80] . The result for the expression of the retarded Green function matrix in the dots is

G r (ε) = 1 D r (ε) g r 1 (ε) g r 1 (ε)Σ r 12 (ε) g r 2 (ε) g r 2 (ε)Σ r 21 (ε) g r 1 (ε) g r 2 (ε) (8)
where

D r (ε) = 1 -g r 1 (ε)Σ r 12 (ε) g r 2 (ε)Σ r 21 (ε), g r i (ε) = g r i (ε)/(1 -Σ r hop,ii (ε)g r i (ε)), g r i (ε) is the retarded Green function in the disconnected dot i defined as g r i (ε) = n∈i g r in (ε) with g r in (ε) = (ε -ε in + i0 + ) -1 , and Σ r ij (ε) is the self-energy given by Σ r ij (ε) = Σ r hop,ij (ε) + δ ji V * ii with Σ r hop,ij (ε) = α=L,R k∈α V * iα g r αk (ε)V jα .
In the multilevel case, it has been assumed that the hopping integrals between the dots, V in,jm , and between the dots and the leads, V in,k , do not depend on the indices n, m and on the state k.

C. Zero-frequency charge susceptibility and dot occupancies

The experimental works carried out on DQDs often focused on establishing the charge stability diagram which gives information on the charge occupancy N i = n∈i d † in d in of each of the two dots 14,27,81,82 . One of the relevant physical quantity to discuss the charge stability diagram is the charge susceptibility which is the linear response in charge Q(t) to the external excitation brought by a time-dependent gate voltage ∆V g (t) applied to the system. Q(t), the charge accumulated on the capacitor plates ensuring the coupling between ∆V g (t) and the dots, is given by 63

Q(t) = (C 0 1 + C 0 2 )∆V g (t) I - i=1,2 α i e∆ N i (t) (9)
where C 0 1 and C 0 2 are the geometrical capacitances of the totally disconnected quantum dots (closed system with V 12 = 0 and V iα = 0, ∀i, α) and I is the identity operator. The last term in Eq. ( 9) comes from the additional electrons in the dot i, denoted by ∆ N i (t), induced by ∆V g (t) when the dots are connected (open system), weighted by the lever-arm coefficient, α i measuring the asymmetry of the capacitive coupling of the voltage generator to the dot i. The external excitation ∆V g (t) introduces the additional source term ∆ H(t) = Q(t)∆V g (t) in the hamiltonian of Eqs. (1)- (3). From the linear response theory [START_REF] Fetter | Quantum Theory of many-particle systems[END_REF] , the expectation value Q(t) , up to the first or-der in ∆V g (t), is given by

Q(t) = (C 0 1 + C 0 2 )∆V g (t) -e 2 ∞ -∞
χ c (t, t )∆V g (t )dt (10) where χ c (t, t ) is the dynamical charge susceptibility given by the Kubo formula:

χ c (t, t ) = i,j=1,2 α i α j χ ij (t, t ), with χ ij (t, t ) = iΘ(t -t ) [∆ N i (t), ∆ N j (t )] (11)
By taking the Fourier transform of Eq. ( 10), one gets

Q(ω) = C(ω)∆V g (ω) where C(ω) = C 0 1 + C 0 2 -e 2 χ c (
ω) is the mesoscopic capacitance 66,67 of the DQD system, with χ c (ω) = i,j=1,2 α i α j χ ij (ω). Thus, in addition to C 0 1 + C 0 2 , C(ω) contains an additional contribution, equal to -e 2 χ c (ω), related to the dynamical charge susceptibility defined as

e 2 χ c (ω) = lim ∆Vg(ω)→0 i=1,2 α i e d ∆ N i (ω) d∆V g (ω) (12)
In the following we will focus on the static charge susceptibility χ c (ω = 0). To get it, we will not make use of Eq. ( 12) but rather of the alternative and more direct expression given by

χ c (ω = 0) = - i,j=1,2 α i α j ∂ N i ∂ε j (13)
where N i is the number of electrons in the dot i at ∆V g (ω) = 0 given by

N i = - i 2π n∈i ∞ -∞ dεG < in,in (ε) (14)
Equation ( 13) can be readily obtained from Eq. ( 12) by noticing that the role of ∆V g (ω = 0) comes down to shift the level energy of the dot j according to ε j → ε j = ε j -α j e∆V g (ω = 0), therefore lim

∆Vg(ω=0)→0 d ∆ N i (ω) d∆V g (ω = 0) = j=1,2 ∂ ∆ N i (ω = 0) ∂ ε j × d ε j d∆V g (ω = 0) = -e j=1,2 α j ∂ N i ∂ε j (15)
Incorporating Eq. ( 15) into Eq. ( 12), one obtains Eq. (13). We set α 1 = α 2 = 1 in the rest of the paper. However, the influence of asymmetric capacitive couplings can be readily studied by using the results we obtain for arbitrary values of α 1 and α 2 .

In this section one has analytically derived all the ingredients needed to characterize the electrical and thermoelectrical properties of a DQD whether it is in serial or parallel geometry. In the next two sections, one successively considers the DQD in series and in parallel for both single-level (SL) and multilevel (ML) dots. The electrical conductance G = dI L /dV , with V = (µ L -µ R )/e the bias voltage between the two leads, and the Seebeck coefficient S = G -1 dI L /d∆T , with ∆T the temperature difference between the left and right leads [START_REF] Crépieux | [END_REF] , are calculated numerically with the help of Eq. ( 4). The zerofrequency charge susceptibility χ c (0) and the total DQD occupancy N = N 1 + N 2 are calculated from Eqs. ( 13) and (14). The calculations are performed in the linear response regime, i.e. in the limit V → 0 and ∆T → 0, but can be readily extended to the nonlinear response regime.

We study the variations of G, S, N , and χ c (0) as a function of gate voltages which act on the positions of energy levels ε 1n and ε 2n of the two dots constituting the DQD system.

III. DISCUSSION FOR A DQD IN SERIES A. SL-DQD in series

We first consider the case of a DQD system in series with a single energy level in each of the two dots (see Fig. 2(a)).

Figure 3 shows the color-scale plots of G, S, N , and χ c (0) as a function of the energies ε 1 and ε 2 , whereas Fig. 4 shows the plots of the same physical quantities as a function of ε 1 along either the first or the second diagonal of equation ε 1 = ε 2 or ε 1 = -ε 2 , respectively. We describe the results obtained in Figs. 3 and4 and then provide an interpretation for them. In a general way, we point out that all the color-scale plots in Fig. 3 have two axes of symmetries which are the first and second diagonals. Figure 3(a) shows that the conductance G is the largest in the central region surrounding the origin point O (ε 1 = 0, ε 2 = 0), with the presence of two peaks along The results obtained for G, S, N and χ c (0) can be interpreted in a simple way by relying on the properties presented in Appendix B. It is explained how, when the DQD system is disconnected from the leads, the hamiltonian H dots describing the system can be diagonalized leading to the eigenenergies E + dots and E - dots , the values of which are given by Eq. (B2). The corresponding antibonding and bonding eigenstates 4,33 |+ and |are defined by Eqs. (B3) and (B4). Consequently the spectral density A 11 (ε) in the dot 1, respectively A 22 (ε) in the dot 2, is a linear combination of delta functions within a multiplicative factor 2π, centered at the values of eigenenergies E + dots and E - dots , corresponding to the spectral densities of eigenstates A + (ε) and A -(ε), with weighting factors equal to |u| 2 and |v| 2 , respectively |v| 2 and |u| 2 in the dot 2, where |u| 2 and |v| 2 are defined by Eqs. (B5) and (B6). The physical meaning of this diagonalization is that the DQD system behaves as an effective single quantum dot with two energy levels at energies E - dots and E + dots . The charge stability diagram of the disconnected DQD system can easily be derived from the latter results. At equilibrium the boundary lines between the domains of different occupations are obtained when any of the two levels of energies E + dots and E - dots is aligned with the chemical potential of the leads, taken to 0 (µ L,R = 0). The equations of the boundary lines B + and B -, hence given by E + dots = 0 and E - dots = 0, are dots is negative while E + dots is positive, hence only the lower energy level at E - dots is occupied. It corresponds either to the domain (1,0) or (0,1) in the charge stability diagram. The line of separation between these latter two domains is along the first diagonal. The results obtained above from the diagonalization of H dots apply to the case of a DQD disconnected from the leads. However it is easy to realize that connecting the DQD system to the leads, would introduce a broadening of the eigenenergy levels described above together with an eventual renormaliza-tion of the eigenenergies values. It would not change the general shape of the charge stability diagram discussed above, but would simply widen the boundary lines separating the different domains. It allows one to explain the charge stability diagram obtained in Fig. 3(c) with a remarkable agreement on the value of the minimal distance observed between the two boundary lines, equal to 2|V 12 |. We now interpret the results obtained for the linear electrical conductance as reported in Figs. 3(a) and 4(a). From Eqs. (A41) and (B3-B6), the linear conductance at zero temperature is proportional to:

ε 1 ε 2 = |V 12 | 2 (see Appendix B
Γ L,11 Γ R,22 |uv| 2 [A + (0) + A -(0)
]. The spectral densities A + (0) and A -(0) are maximal when the point (ε 1 , ε 2 ) corresponding to the energies of the dots falls in one of the boundary lines B + and B -. However in order to get the linear conductance, one has to weight the result for the spectral densities of eigenstate at ε = 0 by a multiplicative coherence factor equal to |uv| 2 . Typically |uv| 2 is the largest along the first diagonal (ε 1 = ε 2 ) where it equals the value 1/4. Moreover u → 1 and v → 0 along the end-part of the boundary lines asymptotic to the ε 2axis, whereas u → 0 and v → 1 along the end-part of the boundary lines asymptotic to the ε 1 -axis. Combining these arguments, the peaks of G in the plane (ε 1 , ε 2 ) arise at the intersection of the boundary lines B + and B -and of the first diagonal as seen in Figs. 3(a) and 4(a). The results obtained for the Seebeck coefficient in Figs. 3(a) and 4(a) can be interpreted in the same way. From Ref. 86, the Seebeck coefficient is proportional to the average energy of the charge carriers: S = -E -µ /k B T which is zero in two types of situations: either when the chemical potential µ (with µ taken to 0 here) falls at the center of one of the peaks of the spectral density of states, that is to say at E + dots and E - dots , or when the chemical potential falls at equal distance from the two peaks. In either case, there are as many carriers with negative energies than carriers with positive energies and hence the average energy of carriers is zero (electron-hole symmetry). This explains why the zeros of S in Fig. 3(b) occur along the boundary lines B -and B + , as well as along the second diagonal, where (E + dots +E - dots )/2 = 0 since ε 1 = -ε 2 there, as seen in Figs. 4(b). The cancellation of S could be used to measure the value of the interdot coupling V 12 which is directly given by the distance between the canceling points. The interpretation of the results obtained for the zero-frequency charge susceptibility follows in the same way: χ(0) determined from Eq. ( 13) is maximal in the plane (ε 1 , ε 2 ) when the spectral densities A + (0) and A -(0) are maximal, i.e. when the point (ε 1 , ε 2 ) falls in one of the boundary lines B + and B -.

B. ML-DQD in series

We examine the case of a ML-DQQ in series. As an example we consider the situation where each of the quantum dots i = 1, 2 has three energy levels of energies ε in = ε i + n∆ε i with the integer index n ∈ [0, 2] and ∆ε i = 1 (see Fig. 2(b)). The results obtained for the color-scale plots of G, S and N , χ c (0) as a function of ε 1 and ε 2 are reported in Figs. 5 and6 respectively in the different interdot coupling regimes. In the weak interdot coupling regime, i.e. for |V 12 | Γ L,11 , Γ R,22 , the electrical conductance G displayed in Fig. 5(a) exhibits peaks centered at the nodes of a square lattice constituted by the vertical lines ε 1n = 0 and horizontal lines ε 2n = 0. Moreover in the continuation of these peaks, one glimpses a slight enhancement of G along the lines of the square lattice. The results for N reported in Fig. 6(a) reveals the charge stability diagram with the presence of 4 × 4 = 16 domains where N changes by plateau. The boundary lines separating the domains coincide with the lines of the square lattice highlighted above. The top-right corner domain corresponds to the completely empty DQD system denoted as (0,0). Whereas one would have expected a completely filled DQD system in the bottom-left corner domain with an occupation (3,3), we point out that the maximal value of N reached there is close to 5, instead of 6, due to the relatively weak value of the dot-lead couplings (Γ L,11 = Γ R,22 = 0.1). In the color-scale plot of χ c (0) shown in Fig. 6(b), χ c (0) is maximal along the same boundary lines as for N , reaching the zero value inside the delimited domains. Finally the color-scale plot of S displayed in Fig. 5(b) shows that S changes of sign several times in the plane (ε 1 , ε 2 ) as previously highlighted in Ref. 87. In the intermediate interdot regime, i.e. for |V 12 | of the order of Γ L,11 , Γ R,22 , the results for the color-scale plots of G displayed in Fig. 5(c) are strongly reminiscent of the results obtained in the case of the single-level DQD in series (see Fig. 3(a)). G is the largest in the central areas surrounding the nodes of the square lattice, with the presence of two peaks on either side of the nodes. These two peaks are located along a line parallel to the first diagonal. When getting farther from these nodes, G gradually decreases along lines forming a characteristic star-shaped pattern, whereas it is zero in the remaining part of the plane. The results for N displayed in Fig. 6(c) still show the presence of 16 domains of different occupancies. The boundary lines between these domains are no longer straight lines but becomes sinuous, the vertices of the square lattice having separated into two triple points in complete agreement with standard results 4,5 . These triple points are at the intersection of the boundaries delimiting three domains of different occupancies. In the color-scale plot of χ c (0) shown in Fig. 6(d), χ c (0) is maximal along the same boundary lines as in Fig. 6(c), reaching the zero value inside the delimited domains. Finally the color-scale plot of S displayed in Fig. 5(d) shows that S exhibits successive changes of sign in the plane (ε 1 , ε 2 ), being of positive sign in the areas delimited by the convex parts of the sinuous boundary line facing top-right, i.e. similar to the area in quadrant I of Fig. 3(b), and of negative sign in its concave parts facing bottom-left, i.e., similar to the area in quadrant III in Fig. 3(b). Besides S is zero inside broad strips surrounding either the second diagonal or the two lines parallel to the second diagonal. Hence S changes sign nine times when sweeping along the first diagonal, instead of five times in Fig. 5(b).

In the strong interdot coupling regime, i.e. for |V 12 | Γ L,11 , Γ R,22 , the square lattice structure visible in the previous figures has disappeared, giving place to an oblique structure in the direction of the second diagonal as can be seen in Figs. 5(e The whole set of these results can be physically interpreted by relying on the properties presented in Appendix B where it is shown how the hamiltonian H dots describing the three-level DQD system can be diagonalized leading to six eigenenergies E λ dots whose values are determined numerically. The charge stability diagram showing N can easily be derived from that. As in the case of the SL-DQD system, the boundary lines between the domains of different occupations are obtained when any of six levels of energies E λ dots is aligned with the chemical potentials of the leads both taken to 0 at equilibrium (µ L,R = 0). The boundary lines of equations E λ dots = 0 are precisely the equienergetic curves represented in Figs. 13(b), 13(d), and 13(f) for different values of the interdot coupling V 12 . One recovers the shapes of the boundary lines obtained above, going from straight lines on a square lattice in the weak interdot coupling regime to sinuous lines in the intermediate and strong interdot coupling regime. As emphasized in Appendix B, as soon as V 12 becomes finite, a level anticrossing effect takes place at the vicinity of the nodes of the square lattice, the vertices of the square lattice splitting into two triple points as observed in Fig. 13. The interpretation of the other results obtained for G, S, and χ c (0) is modeled on the one given above in the single-level case with the presence of some coherence factors which plays a role of extinction for some parts of the equienergetic curves E λ dots = 0.

IV. DISCUSSION FOR A DQD IN PARALLEL A. SL-DQD in parallel

The schematic diagram of the DQD in parallel with a single level in each dot is depicted in Fig. 2(c). The colorscale plots for G, S, N and χ c (0) as a function of ε 1 and ε 2 are shown in Fig. 7. The results are very different from those obtained in the case in series even if in all the plots in Fig. 7, one still glimpses the presence of the boundary lines B + and B -in places almost unchanged compared to the case in series. Strikingly, whereas the color-scale plots of Fig. 3 obtained in the serial cases had two axes of symmetry along the first and the second diagonals, only the axial symmetry with respect to the first diagonal is conserved in Fig. 7 while the one with respect to the second diagonal is lost. The color-scale of G in Fig. 7(a) shows that instead of the two peaks for G in the serial case which occurred at the intersection of the boundary lines B + and B -and of the first diagonal, the positions of the maxima of G now spread all along the boundary line B + in the bottom-left corner, whereas those in the top-right corner are located along the end-parts of the boundary line B -. Moreover one can notice that the conductance ridges thus formed are much broader in the bottom-left corner than in the top-right corner, with higher values reached along the former rather than along the latter ones. We also point out that the amplitude of the conductance to the maximum is about 40% higher than in the serial case. The plots of G along the first and second diagonals displayed in Fig. 8(a) bring a complementary information to that. When sweeping along the first diagonal, G exhibits a single peak at a negative value of ε 1 in perfect agreement with the result displayed in the second diagonal exhibits two peaks located on both sides of the zero value.

As far as the physical interpretation of the results obtained for G, S, N , and χ c (0) is concerned, we would say that even though the basic feature comes from the formation of the boundary lines which occurs when one of the eigenenergies in the system is aligned with the chemical potential of the leads with the coherence factors acting as extinction factors, as explained in Sec. III A, the parallel configuration introduces some noticeable changes with comparison to the serial configuration as the existence of more than one transmission channels. The differences observed in the behaviors of G, S, N , and χ c (0) compared to the case in series results from the interference effects which take place in the presence of two transmission channels as it is the case in the parallel configuration. We give below some simple arguments, developed within the zero temperature and zero interdot coupling limit, which help in elucidating the origin of the formation of the elongated tip mentioned above. In the limit of zero temperature and zero interdot coupling, the expression of N can be derived analytically from Eq. ( 14). One gets for a SL-DQD in series

N series = 1 - 1 π arctan 2ε 1 Γ L,11 - 1 π arctan 2ε 2 Γ R,22 (16)
This result explains why according to the sign of ε 1 and ε 2 , the total occupancy N varies by plateau along which it takes either the value 0, 1 or 2 at most, with a change from one plateau to the other spreading over a width Γ L,11 or Γ R,22 . For a SL-DQD in parallel, one gets when

ε 1 = ε 2 , N parallel = 1 2 - 1 π arctan ε 1 2Γ ( 17
)
where Γ = Γ α,ij for any α = L, R and i, j = 1, 2, i.e., for symmetrical couplings. This expression has to be compared to N series = 1-(2/π) arctan(2ε 1 /Γ) obtained from Eq. ( 16) when ε 1 = ε 2 . In the limit of large negative ε 1 compared to Γ, it leads to N series ≈ 2 and N parallel ≈ 1, explaining the reduction of N by a factor two along the first diagonal (see Fig. 9). Physically, it corresponds to the decoupling of one of the two eigenstates of the DQD, the bonding eigenstate, from the leads.

B. ML-DQD in parallel

We examine the case of a ML-DQQ in parallel schematically represented in Fig. 2

V. CONCLUSION

We have studied a noninteracting DQD system regardless of its geometry, either in series or in parallel, and analyzed its electrical and thermoelectrical properties.

In the case of single-level dots, the expressions for the nonequilibrium Green functions and electrical current are derived exactly. In the case of multilevel dots, an analytical calculation is performed assuming that the hopping integrals between the two dots and between the dots and the leads are independent of energy. The whole set of results apply to any temperatures, bias/gate voltages, and coupling strengths. The numerical results for the linear electrical conductance, the zero-frequency charge susceptibility, the Seebeck coefficient, and the dot occupancy are discussed in the light of previous works.

In particular, the obtained results for G and χ c (0) show that with increasing interdot coupling, the system gradually changes from a regime where the two dots are almost decoupled to a regime where pairs of triple points have separated, until a regime where the two dots merge into a single one, in qualitative agreement with experiments, in particular in the case of a ML-DQD system in series. One observes that for a SL-DQD in series, the Seebeck coefficient undergoes three successive sign changes with increasing dot gates, again in good agreement with experiments. The cancellation of S arises when the average energy of the charge carriers cancels, meaning that the system reaches an electron-hole symmetry situation. A level anticrossing effect resulting from finite interdot coupling manifests itself in the charge stability diagram with boundary lines separating the domains of different dot occupancies occurring when the energy levels of the bonding and antibonding states are aligned with the chemical potentials in the leads. In a striking way, we find a considerable reduction of the total dot occupancy in the case in parallel compared to the case in series, when the energy levels in each of the two dots are equal. We interpret this reduction as an effect of interferences produced by the presence of two transmission electronic paths in the parallel geometry and by the fact that the bonding eigenstate becomes disconnected from the leads. The approach developed in this paper can be directly used to study DQD systems driven out-of-equilibrium by applying either a finite bias voltage or a temperature gradient between the two leads, and/or in the presence of asymmetric couplings. The determination of the noise spectrum and finite-frequency charge susceptibility in the noninteracting DQD system is made in Ref. 88 following the theoretical approach developed in Refs. 89 and 90. A direct and essential extension of this work consists of taking into account the spin degrees of freedom and the Coulomb interactions in the dots.

V in,αk d in (A7) where i = 2 when i = 1 and i = 1 when i = 2. By collecting these contributions together and defining the retarded Green functions in the dots G r in,jm (ε) = d in ; d † jm , one gets the set of coupled equations

G r in,jm (ε) = δ ij δ nm g r in (ε) + g r in (ε) n ∈i V * in,in G r in ,jm (ε) +g r in (ε) α=L,R k∈α V * in,αk G r αk,jm (ε) (A8) G r αk,jm (ε) = g r αk (ε) i =1,2 n ∈i V i n ,αk G r i n ,jm (ε) (A9) where g r in (ε) = (ε -ε in + i0 + ) -1 and g r αk (ε) = (ε - ε αk + i0 + ) -1
are the retarded Green functions in the disconnected dot i and lead α, respectively. By inserting Eq. (A9) into Eq. (A8), one finally obtains a Dyson-like equation

G r in,jm (ε) = δ ij δ nm g r in (ε) + n ∈i g r in (ε)V * in,in G r in ,jm (ε) + α=L,R k∈α i =1,2 n ∈i g r in (ε)V * in,αk g r αk (ε)V i n ,αk G r i n ,jm (ε)(A10)
leading to a set of 4(N ε × M ε ) coupled linear equations. We want to underline that the presence of the term

n ∈i g r in (ε)V * in,in G r in ,jm (
ε) in this equation is directly related to the fact that one has here two coupled dots. In the next two sections, a distinction is made between the SL-DQD case for which an exact solution of Eq. (A10) can be derived and the ML-DQD case for which an approximate solution is given. This latter solution is obtained by making the assumption that the hopping integrals entering in the hamiltonian is independent both of the energy levels ε in in the dots and of the k-state in the leads.

Exact result for a SL-DQD

When each of the two dots contains a single energy level, denoted as ε 1 and ε 2 , the indices n, n and m are absent and then Eq. (A10) reduces to

G r ij (ε) = δ ij g r i (ε) + g r i (ε)V * ii G r ij (ε) + i =1,2 g r i (ε)Σ r hop,ii (ε)G r i j (ε) (A11)
where g r i (ε) = (ε-ε i +i0 + ) -1 is the retarded Green function of the disconnected single-level dot i, and Σ r hop,ii (ε) are the elements of the 2 × 2 hopping self-energy matrix Σ r hop (ε) in the {|1 , |2 }-basis

Σ r hop (ε) = α=L,R k∈α g r αk (ε) |V 1,αk | 2 V * 1,αk V 2,αk V 1,αk V * 2,αk |V 2,αk | 2 (A12)
In a matrix notation, the Dyson equation for G r (ε) writes

G r (ε) = g r (ε) + g r (ε)Σ r (ε)G r (ε) (A13)
where g r (ε) is the retarded Green function matrix of the disconnected dots

g r (ε) = g r 1 (ε) 0 0 g r 2 (ε) (A14)
and Σ r (ε) is the retarded self-energy matrix given by

Σ r (ε) = Σ r hop,11 (ε) Σ r hop,12 (ε) + V * 12 Σ r hop,21 (ε) + V * 21 Σ r hop,22 (ε) (A15)
Equation (A13) can be solved exactly. The calculation presented hereafter applies to the case of a SL-DQD but can be generalized to the case of a ML-DQD within the limit of the approximations made in this paper. The explicit dependences with energy ε of the Green functions and of the self-energies are omitted in the next equations in order to lighten the notations. Starting from Eq. (A11) one writes the four equations of motion followed by the elements of the matrix G r

G r 11 = g r 1 + g r 1 V * 12 G r 21 + g r 1 (Σ r hop,11 G r 11 + Σ r hop,12 G r 21 ) G r 12 = g r 1 V * 12 G r 22 + g r 1 (Σ r hop,11 G r 12 + Σ r hop,12 G r 22 ) G r 21 = g r 2 V * 21 G r 11 + g r 2 (Σ r hop,21 G r 11 + Σ r hop,22 G r 21 ) G r 22 = g r 2 + g r 2 V * 21 G r 12 + g r 2 (Σ r hop,21 G r 12 + Σ r hop,22 G r 22 )
This set of linear equations can be rewritten in the following matrix form M r (ε)G r (ε) = g r (ε), where

M r = 1 -g r 1 Σ r 11 -g r 1 Σ r 12 -g r 2 Σ r 21 1 -g r 2 Σ r 22 (A16) with Σ r ii (ε) = Σ r hop,ii (ε), Σ r ii (ε) = Σ r hop,ii (ε) + V * ii and, G r (ε) = G r 11 (ε) G r 12 (ε) G r 21 (ε) G r 22 (ε) (A17) g r (ε) = g r 1 (ε) 0 0 g r 2 (ε) (A18)
The solution of this matrix equation is given by G r (ε) = (M r (ε)) -1 g r (ε), one obtains the following expression for the 2 × 2 Green function matrix G r (ε)

G r (ε) = 1 D r (ε) g r 1 (ε) g r 1 (ε)Σ r 12 (ε) g r 2 (ε) g r 2 (ε)Σ r 21 (ε) g r 1 (ε) g r 2 (ε) (A19)
where g r i (ε) is defined by

g r i (ε) = g r i (ε)/(1 - Σ r hop,ii (ε)g r i (ε)), D r (ε) is given by D r (ε) = 1 - g r 1 (ε)Σ r 12 (ε) g r 2 (ε)Σ r 21 (ε)
. Equation (A19) gives the exact expression of the retarded Green function G r (ε) in a SL-DQD. It holds as well as for serial as for parallel geometries of the DQD system. We want to underline that for a DQD connected in series, the self-energy Σ r hop (ε) defined in Eq. (A12) becomes a diagonal matrix since the product V * 1,αk V 2,αk is equal to zero whatever the index α is. Therefore, the off-diagonal elements of the total selfenergy matrix Σ r defined in Eq. (A15) reduces to V * 12 and V * 21 in that case.

Generalization to a ML-DQD

In realistic systems, the dots constituting the DQD system contain several energy levels, as for example in Ge/Si heterostructure nanowire-based DQDs 81 or in graphenebased DQDs [START_REF] Li | New Progress on Graphene Research[END_REF] . In that situation one would have to perform a numerical calculation to determine the solutions of Eq. (A10). However, when the hopping integrals V in,im and V in,αk do not depend on the indices n and m and on the state k (and in that case, we use the notations V ii and V iα ), the calculation remains analytical. Within this assumption and by performing a double sum over the n and m indices, Eq. (A10) becomes

G r ij (ε) = δ ij g r i (ε) + g r i (ε)V * ii G r ij (ε) + i =1,2 g r i (ε)Σ r hop,ii (ε)G r i j (ε) (A20)
where

G r ij (ε) = n∈i,m∈j G r in,jm (ε) (A21) g r i (ε) = n∈i g r in (ε) (A22) Σ r hop,ij (ε) = α=L,R k∈α V * iα g r αk (ε)V jα (A23)
In a matrix form, Eq. (A20) reads as

G r (ε) = g r (ε) + g r (ε)Σ r (ε)G r (ε) (A24)
The solutions of Eq. (A24) can be obtained analytically since it is a set of four linear equations. In matrix notation we obtain a 2 × 2 matrix G r (ε), the elements of

which correspond to G r ij (ε), G r (ε) = 1 D r (ε) g r 1 (ε) g r 1 (ε)Σ r 12 (ε) g r 2 (ε) g r 2 (ε)Σ r 21 (ε) g r 1 (ε) g r 2 (ε) (A25) with D r (ε) = 1 -g r 1 (ε)Σ r 12 (ε) g r 2 (ε)Σ r 21 (ε), g r i (ε) = g r i (ε)/(1 -Σ r hop,ii (ε)g r i (ε)) and Σ r ij (ε) = Σ r hop,ij (ε) + δ ji V *
ii . Equation (A25) provides the expression of the retarded Green function G r (ε) in a ML-DQD within the assumption that the hopping integrals are independent of the energy levels. We remark that this result is similar to Eq. (A19) obtained for a SL-DQD, provided that G r (ε), g r i (ε) are changed into G r (ε), g r i (ε). The advanced Green function G a (ε) is obtained straightforwardly by replacing the superscript r by the superscript a in Eq. (A25) with

D a (ε) = 1 -g a 1 (ε)Σ a 12 (ε) g a 2 (ε)Σ a 21 (ε) where g a i (ε) = g a i (ε)/(1 -Σ a hop,ii (ε)g a i (ε)) and Σ a ij (ε) = Σ a hop,ij (ε) + δ ji V ii .
To be able to describe the out-of-equilibrium properties of the DQD such as the electrical current and the electrical conductance, it is necessary to also determine the lesser and greater Green functions G ≶ (ε) for the DQD system as detailed in the next section.

Lesser and greater Green functions

The lesser Green function matrix G < (ε) can be obtained by using the Langreth analytic continuation rules [START_REF] Langreth | [END_REF] on the Dyson equation for the contour ordered Keldysh Green functions obtained from Eq. (A24). One gets

G < (ε) = g < (ε) + g r (ε)Σ r (ε)G < (ε) +g r (ε)Σ < (ε)G a (ε) + g < (ε)Σ a (ε)G a (ε) (A26)
By performing successive iterations on G < (ε) in the r.h.s. of Eq. (A26), one obtains

G < (ε) = G r (ε)Σ < (ε)G a (ε) + 1 + G r (ε)Σ r (ε) g < (ε) 1 + Σ a (ε)G a (ε) (A27)
The second term in the r.h.s of Eq. (A27) vanishes since it can be put in the form

1 + G r (ε)Σ r (ε) g < (ε) 1 + Σ a (ε)G a (ε) = G r (ε)(g r (ε)) -1 g < (ε)(g a (ε)) -1 G a (ε) (A28)
with (g r (ε)) -1 g < (ε)(g a (ε)) -1 = 0, stemming from the fact that g(ε) is the Green function for the disconnected noninteracting DQD system 80 . Therefore, and generalizing it to the greater Green functions G > (ε), one has

G ≶ (ε) = G r (ε)Σ ≶ (ε)G a (ε) (A29)
where the lesser and greater self-energies Σ ≶ (ε) are given by

Σ ≶ (ε) = α=L,R k∈α g ≶ αk (ε) |V 1α | 2 V * 1α V 2α V 1α V * 2α |V 2α | 2 (A30)
Whereas Σ r,a (ε) differs from Σ r,a hop (ε) by the off-diagonal terms V * 12 and V * 21 , Σ ≶ (ε) coincides with Σ ≶ hop (ε); thus one can indifferently use one or the other in any expression where these quantities appear. The result expressed in Eq. (A29) is remarkably simple. It indicates that the information about the inner details of the DQD system is entirely coded in the retarded/advanced Green functions G r,a (ε). We underline that the calculation of G ≶ (ε) is made here for a ML-DQD system. However one can immediately deduce the lesser and greater Green functions G ≶ (ε) for the SL-DQD system by simply changing G r,a (ε) into G r,a (ε) in Eq. (A29). This also applies to the next section.

Explicit expression for the electrical current

The current operator from the lead α is defined as

I α (t) = -ed N α (t)/dt with N α (t) = k∈α c † αk (t) c αk (t).
In the steady state the derivative with respect to the time variable is given by d N α (t)/dt = [ N α (t), H]/i 78 . Thus the average current writes

I α = I α = - e i k∈α [ c † αk (t) c αk (t), H] (A31)
The only term in H leading to a nonvanishing commutator with the product of operators c † αk (t) c αk (t) is H hop . One gets

I α = e k∈α i=1,2 n∈i V iα G < in,αk (t, t) -V * iα G < αk,in (t, t) (A32)
where one has defined the lesser and greater Green functions G < in,αk (t, t ) = i c † αk (t ) d in (t) and G < αk,in (t, t ) = i d † in (t ) c αk (t) . Performing a Fourier transform, one gets

I α = e h k∈α i=1,2 n∈i ∞ -∞ V iα G < in,αk (ε) -V * iα G < αk,in (ε) dε (A33)
In order to calculate the lesser and greater Green functions G < in,αk (ε) and G < αk,in (ε), one applies the Langreth analytic continuation rules [START_REF] Langreth | [END_REF] . From Eq. (A9), one obtains

G < αk,in (ε) = j=1,2 m∈j V jα g r αk (ε)G < jm,in (ε) +g < αk (ε)G a jm,in (ε) (A34) Similarly, G < in,αk (ε) = j=1,2 m∈j V * jα G < in,jm (ε)g a αk (ε) +G r in,jm (ε)g < αk (ε) (A35)
By inserting these expressions into Eq. (A32), one gets

I α = e h k∈α i=1,2 j=1,2 ∞ -∞ V iα V * jα × g a αk (ε) -g r αk (ε) G < ij (ε) +g < αk (ε) G r ij (ε) -G a ij (ε) dε (A36)
where G r ij (ε) is the Green function summed over the indices n and m as defined in Eq. (A21). By using the general relationship G r (ε) -G a (ε) = G > (ε) -G < (ε) which holds for any Green function G, one obtains

I α = e h k∈α i=1,2 j=1,2 ∞ -∞ V iα V * jα × g < αk (ε)G > ij (ε) -g > αk (ε)G < ij (ε) dε (A37)
Physically, this expression is interpreted as follows: the first contribution in I α represents the current flowing from the α lead to the DQD since it is the product of the out-tunneling rate of the occupied state in the α lead, k∈α V iα V * jα g < αk which corresponds to the self-energy, and of the number of unoccupied states in the DQD, G > ij (ε), whereas the second contribution with the minus sign corresponds to the current flowing from the DQD to the lead α. Equation (A37) can be written thanks to Eq. (A29) under the following form

I α = e h ∞ -∞ Tr Σ < α (ε)G r (ε)Σ > (ε)G a (ε) -Σ > α (ε)G r (ε)Σ < (ε)G a (ε) dε (A38)
where the matrix elements of the self-energy

Σ ≶ α (ε) are defined as Σ ≶ α,ij (ε) = k∈α V * iα g ≶ αk (
ε)V jα and where Tr[ ] denotes the trace of the matrix. In the limit of wide flat band for elecrons in the leads and energy-independent hopping integrals, one has

Σ < α (ε) = if α (ε)Γ α (A39) Σ > α (ε) = -i(1 -f α (ε))Γ α (A40)
where the elements Γ α,ij of the matrix Γ α are defined as Γ α,ij = 2πρ α V * iα V jα with ρ α , the density of states in the lead α. We also have Σ r,a α (ε) = ∓iΓ α /2. By inserting Eqs. (A39) and (A40) into Eq. (A38), one obtains

I α = e h ∞ -∞ T αα (ε) f α (ε) -f α (ε) dε (A41)
which corresponds to Eq. ( 4), where T αα (ε) is the transmission coefficient equal to

T αα (ε) = Tr Γ α G r (ε) Γ α G a (ε) (A42)
with α = R for α = L and α = L for α = R, and where Γ α is the dot-lead coupling matrix defined as For a SL-DQD, the hamiltonian H dots writes as a 2 × 2 matrix in the basis {|1 , |2 } of the states in the two dots 1 and 2

Γ α = 2πρ α |V 1α | 2 V * 1α V 2α V 1α V * 2α |V 2α | 2 (A43)
H dots = ε 1 V * 12 V * 21 ε 2 (B1)
It can be diagonalized leading to the following eigenenergies E + dots and E - dots

E ± dots = ε 1 + ε 2 ± (ε 1 -ε 2 ) 2 + 4|V 12 | 2 2 (B2)
and eigenvectors |+ and |-which correspond to the antibonding and bonding eigenstates of the SL-DQD

|+ = u * |1 + v * |2 (B3) |-= -v|1 + u|2 (B4)
with

|u| 2 = 1 2 1 + ε 1 -ε 2 E + dots -E - dots (B5) |v| 2 = 1 2 1 - ε 1 -ε 2 E + dots -E - dots (B6)
and therefore |uv| 2 = |V 12 | 2 /(E + dots -E - dots ) 2 . It gives rise to a level anticrossing effect as soon as the interdot coupling V 12 is finite. The anticrossing of the two levels ε 1 and ε 2 occurs in the vicinity of the first diagonal as shown in Fig. 12. Indeed, from Eq. (B2), the difference between the two eigenenergies reads as ∆E = E + dots -E - dots = (ε 1 -ε 2 ) 2 + 4|V 12 | 2 , which is minimal along the first diagonal (i.e., when ε 1 = ε 2 ) equaling ∆E = 2|V 12 | then. Resulting from Eqs. (B3)-(B6), the spectral density A 11 (ε) in dot 1, respectively, A 22 (ε) in dot 2, is a linear combination of Dirac delta functions within a multiplicative factor 2π, centered at the values of eigenenergies E + dots and E - dots , with weighting factors equal to |u| 2 and |v| 2 , respectively, |v| 2 and |u| 2 in dot 2. It has to be noted that as soon as the interdot coupling V 12 becomes finite, a mixed spectral density A 12 (ε) arises, resulting from interdot transitions. The charge stability diagram of the system can easily be derived from the latter results. At equilibrium the boundary lines between the domains of different occupations are obtained when any of the two levels of energies, E + dots and E - dots is aligned with the chemical potential of the leads µ L = µ R = µ. The equations of the boundary lines B + and B -, hence given by E + dots = µ and E - dots = µ, are (ε 1 -µ)(ε 2 -µ) = V 2 12 . They correspond to two branches of an hyperbol in the plane (ε 1 , ε 2 ).The distance between the two branches B + and B -is minimal along the first diagonal taking the value of 2|V 12 |. The boundary lines are drawn in Figs. 13(a), 13(c), and 13(e) at µ = 0 for different values of the interdot coupling V 12 . One can check that the minimal distance between the two boundary lines B -and B + increases with increasing V 12 . When V 12 = 0, the eigenenergies are simply equal to ε 1 , ε 1 + ∆ε 1 , ε 1 + 2∆ε 1 , ε 2 , ε 2 + ∆ε 2 and ε 2 + 2∆ε 2 , and the equienergetic curves at E λ dots = 0 are the three horizontal and three vertical lines of a square lattice as found in Fig. 13(b). As soon as the interdot coupling V 12 gets finite, a level anticrossing effect takes place in the vicinity of the nodes of the square lattice, as shown in Figs. 13(d) and 13(f). At V 12 = 0.1, the equienergetic curves become sinuous, as a result of this level anticrossing effect. At V 12 = 0.5, the distance between two adjoining equienergetic curves increases, leaving room for wide interstitial areas in the direction parallel to the second diagonal. These various elements brought by the above discussion are crucial to physically interpret the results obtained for the conductance, Seebeck coefficient, total dot occupancy, and zero-frequency charge susceptibility in the DQD system as discussed in Secs. III A and III B.

FIG. 1 .

 1 FIG. 1. Schematic representation of the DQD coupled to L and R leads in the serial (a) and parallel (b) geometries.

FIG. 2 .

 2 FIG. 2. The four studied geometries: (a) SL-DQD in series, (b) ML-DQD in series, (c) SL-DQD in parallel, and (d) ML-DQD in parallel. We take ε1n = ε1 + n∆ε1 and ε2n = ε2 + n∆ε2 with the integer index n ∈ [0, 2] in panels (b) and (d) since one considers three energy levels in each of the two dots in the ML-DQD case. The dotted black lines symbolize the couplings between the various parts of the system.

FIG. 3 .FIG. 4 .

 34 FIG. 3. Color-scale plots of (a) the linear conductance G, (b) the Seebeck coefficient S, (c) the total dot occupancy N , and (d) the zero-frequency charge susceptibility χc(0) for a SL-DQD connected in series as a function of ε1 and ε2 for µL,R = 0, kBT = 0.01, V12 = 0.1, and ΓL,11 = ΓR,22 = 0.1. In panel (a) the black line marks the first diagonal of equation ε1 = ε2, whereas the red line shows the second diagonal of equation ε1 = -ε2 along which the plots of Fig. 4 are drawn. In panel (c) the domain with occupancies N1 and N2 of the dots 1 and 2 are indicated under the form ( N1 , N2 ), and the dashed and dotted black arcs represent the boundary lines B+ and B-between domains with different occupancies. Panel (d) shows the four quadrants: I (top right), II (top left), III (bottom left) and IV (bottom right).

  ). They correspond to the two branches of an hyperbole as shown in Figs.13(a), (c), and (e) for different values of the interdot coupling V 12 . The distance between the two branches is minimal along the first diagonal, taking the value of 2|V 12 |. They correspond to the two boundary lines found in quadrants I and III in Fig. 3(c): the curve B -of equation E - dots = 0 corresponds to the boundary line in quadrant I, whereas the curve B + of equation E + dots = 0 corresponds to the boundary line in quadrant III. Inside the top-right domain delimited by the curve B -, both E - dots and E + dots are positive and the two corresponding energy levels are empty. It gives rise to the domain (0,0) in the charge stability diagram (by making use of the results on the spectral densities A ii (ε) mentioned above and by noticing that |u| 2 +|v| 2 = 1). Inside the bottom-left domain delimited by the curve B + , both E - dots and E + dots are negative and the corresponding two energy levels are occupied. It gives rise to the domain (1,1) in the charge stability diagram. Inside the area located between the two boundary lines B + and B -, centered around the second diagonal, E -

  ) and 5(f), and Figs. 6(e) and 6(f) for G, S, N and χ c (0). This means that the two quantum dots have merged into one single quantum dot of occupation N 1 + N 2 . In Fig.6(e), one sees that the change of N by plateau observed in the weak and intermediate regimes, is replaced by a smooth variation. The orders of magnitude obtained for N and χ c (0) are ten times smaller in the case of strong coupling regime compared to the weak and intermediate coupling regimes, the observed reduction being of the same magnitude as the reduction of Γ L,11 and Γ R,22 values. The results for S shown in Fig.5(f) are in qualitative agreement with the experimental ones obtained in GaAs/AlGaAs heterostructures displayed in Fig.2of Ref.58.

FIG. 5 .

 5 FIG. 5. Color-scale plots of G and S for a ML-DQD connected in series as a function of ε1 and ε2 at µL,R = 0, kBT = 0.01 for: (a), (b) weak interdot coupling (V12 = 0.01, ΓL,11 = ΓR,22 = 0.1); (c), (d) intermediate interdot coupling (V12 = ΓL,11 = ΓR,22 = 0.1); and (e), (f) strong interdot coupling (V12 = 0.5, ΓL,11 = ΓR,22 = 0.01).

FIG. 6 .

 6 FIG.6. Color-scale plots of N and χc(0) for a ML-DQD connected in series as a function of ε1 and ε2. The parameters are the same as in Fig.5.

Fig. 2 of

 2 FIG. 7. Color-scale plots of (a) the linear electrical conductance G, (b) the Seebeck coefficient S, (c) the total dot occupancy N , and (d) the zero-frequency charge susceptibility χc(0) for a SL-DQD connected in parallel as a function of ε1 and ε2 for µL,R = 0, kBT = 0.01, V12 = 0.1, and Γα,ij = 0.1 for both α = L, R and i, j = 1, 2.

  (d) taking as an example the case of three energy levels of energies as in the case of the ML-DQD in series. The results obtained for the colorscale plots of G, S and N , χ c (0) as a function of ε 1 and ε 2 are reported in Figs.10 and 11, respectively, in the different interdot coupling regimes. In the weak interdot coupling regime, i.e., for |V 12 | Γ α,ij , the color-scale plot for G reported in Fig.10(a) is strongly modified compared to the configuration in series, showing conductance ridges along the lines of a slightly distorted square lattice. The color-scale plot of χ c (0) shown in Fig. 11(b) looks like that of G, with maxima along the lines of a lattice, with the presence of additional pairs of localized peaks around the nodes of the lattice. The results for N reported in Fig. 11(a) reveal the charge stability diagram with the presence of 4×4 = 16 domains as for the configuration in series. However one notices an important difference which is provided by the presence of elongated tips along the first diagonal as well as along the secondary first diagonals, with a dip of N along them, similarly to what is observed in the case of the SL-DQD in parallel. Finally the color-scale plot of S displayed in Fig. 10(b) shows that S changes sign several times in the plane (ε 1 , ε 2 ). In the intermediate interdot regime, i.e., for |V 12 | of the order of Γ α,ij , the results for the different color-scale plots of G, S, N and χ c (0) are strongly reminiscent of the results reported in Fig. 7 obtained in the case of the SL-DQD in parallel with the same choice of parameters. The reported pattern corresponds to the duplication of the pattern observed in the single-level case, in each cell of a square lattice. Here again the presence of the elongated tips are clearly visible in the charge stability diagram revealed by the color-scale plot of N displayed in Fig. 11(d). In the strong interdot coupling regime, i.e., for |V 12 | Γ α,ij , the square lattice structure visible in the previous figures has disappeared, giving place to an oblique structure in the direction of the second diagonal as can be seen in Figs. 10(e) and 10(f), and 11(e) and 11(f) for G, S, N , and χ c (0) respectively. In addition, the extremities of the lines of maxima for G, S, and χ c (0) disappear in the top-right part of the plane (ε 1 , ε 2 ) compared to the case of ML-DQD in series. The color-scale plot of N shows the presence of a predominant tip along the first diagonal.

FIG. 10 .

 10 FIG. 10. Color-scale plots of G and S for a ML-DQD connected in parallel as a function of ε1 and ε2 at µL,R = 0, kBT = 0.01 for: (a), (b) weak interdot coupling (V12 = 0.01 and Γα,ij = 0.1); (c), (d) intermediate interdot coupling (V12 = Γα,ij = 0.1); and (e), (f) strong interdot coupling (V12 = 0.5, Γα,ij = 0.01).

FIG. 11 .

 11 FIG. 11. Color-scale plots of N and χc(0) for a ML-DQD connected in parallel as a function of ε1 and ε2. The parameters are the same as in Fig. 10.

FIG. 12 .

 12 FIG.12. Eigenenergies E + dots (red curve) and E - dots (black curve) at V12 = 0.1 as a function of the angle θ with ε1 = cos(θ) and ε2 = sin(θ). Level anticrossing effect arises at θ = π/4 or θ = 5π/4, i.e. at ε1 = ε2, with a distance between the red and black curves equal to ∆E = 2|V12| = 0.2.

2 .

 2 ML-DQD For the three energy level ML-DQD considered in Sec. III B, the hamiltonian H dots writes as a 6 × 6 matrix in the basis of the states {|1n , |2m } with integer indices n, m ∈ [0, 2] in the two dots given by H dots = where ε in = ε i + n∆ε i and V 12 is the interdot coupling of equal value regardless of the levels considered in each dot. H dots can be diagonalized leading to six eigenenergies E λ dots , with the integer index λ ∈ [1, 6], whose values can be numerically calculated. It is useful for the discussion in Sec. III B to determine the equienergetic curves at equilibrium: E λ dots = µ. The results obtained numerically are displayed in Figs. 13(b), 13(d), and 13(f) at µ = 0 for different values of the interdot coupling V 12 .

FIG. 13 .

 13 FIG. 13. Equienergetic curves of equations E ± dots = µL,R (left column) and E λ dots = µL,R with λ ∈ [1, 6] (right column) at µL,R = 0 in the plane (ε1, ε2) for (a)-(b) V12 = 0, (c)-(d) V12 = 0.1 and (e)-(f) V12 = 0.5. The purple lines are obtained for a SL-DQD and the black lines for a ML-DQD with three levels of energy in each dot.

Acknowledgments -The authors would like to thank H. Baranger, S. Barraud, A. Crippa, A.-M. Daré, S. De Franceschi, V. Kashcheyevs, R. Maurand, C. Mora, Y.-M. Niquet, S. Sahoo, M. Sanquer and R. Whitney for valuable discussions. They acknowledge the CEA-Programme Transverse Nanosciences and CEA-Eurotalents for financial support.

Appendix A: Electrical current in a DQD

Equations of motion

The Green functions in the dots are obtained by using the equation of motion approach. We adopt the Zubarev notation 91 along which the retarded Green function in energy G r A, B (ε) associated with operators A and B is denoted by A; B , with G r A, B (ε) the Fourier transform of the retarded Green function in time, G r A, B (t, t ) = -Θ(t -t ) { A, B} , Θ being the Heaviside step function. By using this notation, the equation of motion for A; B writes

where [ A, H] is the commutator between the operator A and the hamiltonian H of the DQD system, and { A, B} is the anticommutator between the operators A and B. Applying Eq. (A1) to the various operators in the dots and the leads, one gets

The calculation of the nonvanishing commutators between the operators d in and c αk and the various terms of H leads to

Appendix B: H dots eigenvalues

In this appendix, we determine the eigenenergies and eigenvectors of the hamiltonian H dots of Eq. ( 1) describing the DQD disconnected from the leads firstly for a SL-DQD and secondly for a ML-DQD with three levels of energy.