{'fraction_non_alphanumeric': 0.0450286458832321, 'fraction_numerical': 0.015827089999808388, 'mean_word_length': 4.014893821466321, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Level anticrossing effect in single-level or multi-level double quantum dots: Conductance, zero-frequency charge susceptibility and Seebeck coefficient

A theory is developed to model a double quantum dot in a unified way whatever the geometry of the two dots is, either in series or in parallel. In the case of a single-level double quantum dot, the problem is exactly solvable whereas for a multi-level double quantum dot, an analytical solution is obtained in the limit of energy-independent hopping integrals. After deriving the expressions for the non-equilibrium Green functions, we study the dependences of the conductance, zero-frequency charge susceptibility and Seebeck coefficient on the gate voltages applied to the dots, allowing us to derive the charge stability diagram revealing a wide variety of behaviors for the system. The findings are in agreement with the experimental observations notably with the occurrence of successive sign changes of the Seebeck coefficient when varying the gate voltages. We interpret the results in terms of the bonding and antibonding states produced by the level anticrossing effect which occurs in the presence of a finite interdot coupling. We show that at equilibrium the boundary lines between the domains with different dot occupancies in the charge stability diagram, take place when the bonding and antibonding state levels are aligned with the chemical potentials in the leads. Finally the total dot occupancy is found to be considerably reduced in the case in parallel compared with the case in series, when the level energies in each dot are equal. We interpret this dip as a direct manifestation of the interference effects occurring in the presence of the two electronic transmission paths provided by each dot.


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 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 to 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] . The serial-coupling 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 [29][30][31] or other specific effects [32][33][34][35][36][37] . In this article, both geometries of serial-coupled and parallel-coupled DQDs are considered.
From the theoretical side, the electrical transport properties in DQD have been widely studied and the approaches to model them are mainly based on the use of Master equations for the density matrix [38][39][40][41][42][43] or on the use of Keldysh Green function methods [44][45][46][47][48][49] which becomes es-sential when one wants to treat out-of-equilibrium situations. The overall evolution of the linear conductance as a function of gate voltages is well understood within a classical theory 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-6 : (i) at weak interdot coupling, conductance peaks are observed at the nodes of a square lattice, (ii) at intermediate interdot coupling, pairs of triple-points arise at the boundaries between the regions with different electron occupancy in the dots to form a honeycomb lattice in the plane (ε 1 , ε 2 ), where ε 1 and ε 2 are the energy levels of the two dots, 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 number of electrons in the dots 1 and 2 respectively.
Theoretical works devoted to the thermoelectrical properties of DQDs have only started in the last ten years. Various approaches have been used such as Master equations [50][51][52] , non-crossing approximation 53 , numerical renormalization group 54 , non-equilibrium Green function approach 55,56 or linear-response theory [57][58][59] . At the same time, experimental studies on thermopower, i.e. Seebeck coefficient, are currently growing up [60][61][62][63][64] , so that a continuous effort must be provided simultaneously on the theoretical side.
In this paper, we develop theoretical tools to describe both the electrical and thermoelectrical properties of single-level (SL) and multi-level (ML) DQDs connected to two reservoirs (leads) of electrons. The determination of the non-equilibrium Green functions in the DQD system allows us to derive the analytical expressions for the electrical current and the dot occupancy. The results are valid at any temperature, lead-dot and interdot couplings, bias and gate voltages. The assumptions made in this work are the following: (i) we neglect Coulomb interactions, (ii) we consider the wide flat band limit for electrons in the leads, and (iii) we assume in the ML-DQD case that the various couplings are independent on the values of the energy levels in the dots. The obtained results apply to both parallel and series geometries of the DQD. We numerically calculate the conductance, the dot occupancies, the zero-frequency charge susceptibility and the Seebeck coefficient in various regimes going from weak to strong interdot coupling. The results for the evolution of these quantities as a function of gate voltages are in qualitative agreement with experiments. This paper is structured as follows: in Sec. II, the hamiltonian used to describe the DQD system is given and the calculation of the Green functions for both SL-DQD and ML-DQD are detailed. The determination of the electrical current and of the dot occupancy is respectively made in Secs. III and IV. The numerical results for the conductance, the Seebeck coefficient, the dot occupancy and the zero-frequency charge susceptibility are presented for a DQD connected in series in Sec. V, and for a DQD connected in parallel in Sec. VI. We conclude in Sec. VII.

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 respectively with n ∈ [0, N ε −1] and m ∈ [0, M ε −1]. The values of ε 1n and ε 2m can be tuned by varying the nearby gate voltages. The two metallic left (L) and right (R) leads are characterized by their chemical potentials µ L , µ R and temperatures T L , T R respectively. In the non-interacting case that we consider, the hamiltonian of this DQD system writes: H = H dots + H leads + H hop , where H dots and H leads are the hamiltonian of the disconnected DQD and that of the disconnected leads respectively, and H hop is the hopping hamiltonian between the dots and the leads 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 character of the expression we introduce for the hamiltonian. It allows one to describe all the possible assemblies of two dots in an 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 non-interacting DQD system connected to non-magnetic leads.

B. Equations of motion
The Green functions in the dots are obtained by using the equation of motion approach. We adopt the Zubarev notation 65 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 The calculation of the non-vanishing commutators between the operators d in and c αk and the various terms of H leads to where i = 2 when i = 1 and i = 1 when i = 2. By collecting these contributions together and defining the retarded Green function G r in,jm (ε) = d in ; d † jm , one gets the set of coupled equations 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. (12) into Eq. (11), one finally obtains a Dyson-like equation (13) 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. (13) 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.

C. 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. (13) reduces to 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 In a matrix notation, the Dyson equation for G r (ε) writes G r (ε) = g r (ε) + g r (ε)Σ r (ε)G r (ε) (16) where g r (ε) is the retarded Green function matrix of the disconnected dots g r (ε) = g r 1 (ε) 0 0 g r 2 (ε) (17) and Σ r (ε) is the retarded self-energy matrix given by Eq. (16) can be solved exactly. The details of the calculation are given in Appendix A. We obtain the following expression for the 2 × 2 Green function matrix G r (ε) 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 (ε). Eq. (19) 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. (15) 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. (18) reduces to V * 12 and V * 21 in that case.

D. 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 66 or in graphenebased DQDs 67 . In that situation one would have to perform a numerical calculation to determine the solutions of Eq. (13). 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. (13) becomes where In a matrix form, Eq. (20) reads as The solutions of Eq. (24) 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 (ε), ii . Eq. (25) provides the expression of the retarded Green function G r (ε) in a ML-DQD within the assumption that the hopping integrals in independent of the energy levels. We remark that this result is similar to Eq. (19) 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. (25) 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 current and the conductance, it is necessary to also determine the lesser and greater Green functions G ≶ (ε) for the DQD system as detailed in the next section.

E. Lesser and greater Green functions
The lesser Green function matrix G < (ε) can be obtained by using the Langreth analytic continuation rules 68 on the Dyson equation for the contour ordered Keldysh Green functions obtained from Eq. (24). One gets By performing successive iterations on G < (ε) in the r.h.s. of Eq. (26), one obtains The second term in the r.h.s of Eq. (27) vanishes since it can be put in the form = G r (ε)(g r (ε)) −1 g < (ε)(g a (ε)) −1 G a (ε) (28) with (g r (ε)) −1 g < (ε)(g a (ε)) −1 = 0, stemming from the fact that g(ε) is the Green function for the disconnected non-interacting DQD system 69 . Therefore, and generalizing it to the greater Green functions G > (ε), one has where the lesser and greater self-energies Σ ≶ (ε) are given by 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. (29) 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. (29). This also applies to the next section.

III. ELECTRICAL CURRENT
In this section we derive the expression for the electrical current assuming that the non-interacting DQD is in a steady state. For that we start from the current operator from the lead α defined as In the steady state the derivative with respect to the time variable is given by 70 , H]/i . Thus the average current writes The only term in H leading to a non-vanishing commutator with the product of operators c † αk (t) c αk (t) is H hop . One gets where one has defined the lesser and greater Green func- Performing a Fourier transform, one gets In order to calculate the lesser and greater Green functions G < in,αk (ε) and G < αk,in (ε), one applies the Langreth analytic continuation rules 68 . From Eq. (12), one obtains Similarly, By inserting these expressions into Eq. (32), one gets where G r ij (ε) is the Green function summed over the indices n and m as defined in Eq. (21). By using the general relationship G r (ε) − G a (ε) = G > (ε) − G < (ε) which holds for any Green function G, one obtains 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 α. Eq. (37) can be written thanks to Eq. (29) under the following form where the matrix elements of the self-energy Tr[ ] denotes the trace in the basis of the states |1n and |2m . The calculation leads to where the elements Γ α,ij of the matrix Γ α are defined as Γ α,ij = 2πρ α V * iα V jα . We also have Σ r,a α (ε) = ∓iΓ α /2. We underline that the results above are obtained in the wide flat band limit neglecting the energy dependency of the density of states ρ α of the lead α. By inserting Eqs. (39) and (40) into Eq. (38), one obtains where α = R for α = L and α = L for α = R. This result explicitly shows that the current from the left lead and the current from the right lead are opposite in sign: meaning that the electrical current is a conserved quantity as expected in the steady state. Eq. (41) corresponds to the Landauer formula for the electrical current with a transmission coefficient equal to where Γ α is the dot-lead coupling matrix defined as This latter matrix is diagonal for a DQD in series but non-diagonal for a DQD in parallel.

IV. 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,66,71 . 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 72 where C 0 1 and C 0 2 are the 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. (44) 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 73 , the expectation value Q(t) , up to the first order in ∆V g (t), is given by where χ c (t, t ) is the dynamical charge susceptibility given by the Kubo formula: By taking the Fourier transform of Eq. (45), one gets contains an additional contribution, equal to −e 2 χ c (ω), related to the dynamical charge susceptibility defined as In the following we will focus on the static charge susceptibility χ c (ω = 0). To get it, we will not make use of Eq. (47) but rather of the alternative and more direct expression given by where N i is the number of electrons in the dot i at ∆V g (ω) = 0 given by Eq. (48) can be readily obtained from Eq. (47) 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 Incorporating Eq. (50) into Eq. (47), one obtains Eq. (48). 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 Secs. II, III and IV, one has derived analytically 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 SL and ML dots. The 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 76 , are calculated numerically with the help of Eq. (41). The zero-frequency charge susceptibility χ c (0) and the total DQD occupancy N = N 1 + N 2 are calculated from Eqs. (48) and (49). 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 non-linear 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.

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)). Fig. 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 and 4 and then provide an interpretation for them. In a general way, we point out that all the color-plots in Fig. 3 have two axes of symmetries which are the first and second diagonals. Fig. 3(a) shows that the conductance G is the  Fig. 4(a)). When getting farther from the origin O, G gradually decreases forming a star-shaped pattern in the plane (ε 1 , ε 2 ) as displayed in blue color in Fig. 3(a), until reaching the zero value in the remaining parts of the plane. The color-plot of the Seebeck coefficient displayed in Fig. 3(b) shows that S is zero (green color) along the boundary lines B + and B − located in quadrants III and I, as well as inside a band located on either side of the second diagonal. It takes positive values (red color) inside the top-right domain delimited by the boundary line B − , and vice-versa negative values (violet color) inside the respective bottom-left domain delimited by B + . In the intermediate area located between the two boundary lines, S exhibits pockets of local maxima and minima located in the vicinity of B + and B − . In the plot of S as a function of ε 1 along the first and second diagonal displayed in Fig. 4(b), one can see that S is zero all along the second diagonal, while it changes of sign three times when sweeping along the first diagonal, once at ε 1 = 0 and the two other times at the positions of the maxima of G observed in Fig. 4(a). The behavior for S which we obtain with these three observed changes of sign, are in agreement with the experimental results obtained in GaAs/AlGaAs heterostructures which are reported in Fig. 3 of Ref. 60 and Fig. 4 in Ref. 63. We underline that in the experimental works V th = −S∆T is plotted instead of S in the results presented here. Fig. 3(c) shows the evolution of N as a function of ε 1 and ε 2 . It reveals the charge stability diagram with the presence of four domains denoted as ( N 1 , N 2 ) inside which N 1 and N 2 take values close to (0,0), (1,0), (0,1) or (1,1). As can be seen, the (0,0) and (1,1) domains are shrunk within the I and III quadrants, compared to the uncoupled DQD case (at V 12 = 0) for which the diagram would have shown a tiling on a square lattice set on the four quadrants (see Figs. 13(a), (c) and (e) in Appendix B). Note that the boundary lines B − and B + delimiting the (0,0) and (1,1) domains have the shape of arcs and that N takes a constant value along the second diagonal (see red curve in Fig. 4(c)). The evolution of the zero-frequency charge susceptibility χ c (0) displayed in Figs. 3(d) and 4(d) follows the same trend. It shows the presence of two lines of maxima, located precisely on the boundary lines B + and B − highlighted in Fig. 3(c). As for G, the origin O behaves as a saddle point, χ c (0) being maximal when sweeping along the second diagonal direction, while it is a local minimum along the first diagonal one.
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 |+ 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 ε 1 ε 2 = |V 12 | 2 (see Appendix B). 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 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 − 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 renormalization 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 conductance as reported in Figs. 3(a) and 4(a). From Eqs. (41) and (B3-B6), the linear conductance at zero temperature is proportional to: Γ 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 ε 2 -axis, 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. 77, 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 (particle-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 interpretation of the results obtained for the zero-frequency charge susceptibility follows in the same way: χ(0) determined from Eq. (48) 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-plots of G, S and N , χ c (0) as a function of ε 1 and ε 2 are reported in Figs. 5 and 6 respectively in the different interdot coupling regimes.
In the weak interdot coupling regime, i.e. for |V 12 | Γ L,11 , Γ R,22 , the 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 topright 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 bottomleft 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 dotlead couplings (Γ L,11 = Γ R,22 = 0.1). In the color-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-plot of S displayed in Fig. 5(b) shows that S changes of sign several times in the plane (ε 1 , ε 2 ). In the intermediate interdot regime, i.e. for |V 12 | of the order of Γ L,11 , Γ R, 22 , the results for the color-plots of G displayed in Fig. 5(c) is 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 starshaped pattern, whereas it is zero in the remaining part of the plane. The results for N displayed in Fig. 6(c) still shows 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. These triple points are at the intersection of the boundaries delimiting three domains of different occupancies. In the color-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-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)-(f) and Figs. 6(e)-(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. 2 of Ref. 60. 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 diagonal-ized leading to six eigenergies 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), (d) and (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.

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 is much broader in the bottom-left corner than the one in the top-right corner, with higher values reached along the former 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 , whereas the peak of G along the second diagonal is centered at the zero value. The results obtained for S displayed in Fig. 7(b) shows that S still takes positive values inside the top-right domain and negative values inside the bottom-left domain as it was the case for the single-level DQD case, with a strong reduction in the order of magnitude of the amplitude compared to the serial configuration results shown in Fig. 3(b). In the intermediate area located between these two domains, S exhibits a series of minima of negative sign inside a pair of two triangular pockets along the boundary line delimiting the top-right domain as in Fig. 3(b). Nevertheless we point out three major differences: (i) the base of these triangles is a straight line parallel to the second diagonal shifted to the top-right corner from the second diagonal by a distance equal to |V 12 | and one has S = 0 on this base; (ii) the two triangles are disjoint at the center showing a gap around the first diagonal in which S takes a positive value close to zero; and (iii) the negative values reached inside the triangles are about twice larger than inside the bottom-left corner. Moreover one still observes the presence a series of maxima of positive sign inside a pair of two triangular pockets as in Fig. 3(b), the difference in the parallel case is that these two triangular pockets are contiguous to the previous ones. Besides they share with the negative sign triangles the same three peculiarities i.e. the area located between the latter positive sign triangles and the boundary line surrounding the bottom-left domain, exhibiting a gap around the first diagonal. The plots of S along the first and second diagonals displayed in Fig. 8(b) completes these information. Along the first diagonal, S goes from negative to positive values with increasing ε 1 , showing one change of sign instead of three in the serial configuration. Along the second diagonal, S is no longer zero since the electron-hole symmetry holding in the serial case is now lost, instead S keeps a positive value all along the second diagonal with a marked minimum around ε 1 = 0. The electron-hole symmetry is however restored with S = 0 on a line which corresponds to the base of the triangles discussed there before (see the green line parallel to the second diagonal in Fig. 7(b)). The color-plot of N displayed in Fig. 7(c) still shows the presence of the boundary lines in top-right and bottomleft corners as in Fig. 3(c). However one can notice important differences in comparison to the results obtained for the case in series: (i) the value of N is strongly reduced inside the domain in the bottom-left corner, reaching the value 0.8 instead of 1.6 ; (ii) an elongated tip is formed along the first diagonal extending from the domain in top-right corner to the other side where it cuts the bottom-left corner domain in half. Strikingly the value of N is strongly reduced along this tip going from 0 to 0.8 when sweeping along the first diagonal as shown in Fig. 8(c). Moreover the plot of N along the second diagonal displayed in Fig. 8(c) too shows a marked minimum at ε 1 = 0 at the crossing with the elongated tip previously reported. Finally the evolution of χ c (0) displayed in Fig. 7(d) follows the same trend. It shows lines of maxima along the same boundary lines highlighted in Fig. 7(a) but the maxima are much more pronounced along the boundary lines in the top-right corner in comparison to the ones in the bottom-left corner. Besides these two lines of maxima are cut into half at the crossing with the first diagonal with the opening of a gap around it. As shown in Fig. 8(d), when sweeping along the first diagonal, χ c (0) exhibits a peak at a negative sign value of ε 1 , whereas χ c (0) along 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. V 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. (49). One gets for a SL-DQD in series 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 , 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. (51) 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(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 colorplots 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-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-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) reveals 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-plot of S displayed in Fig. 10(b) shows that S changes of 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-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-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)-(f) and 11(e)-(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-plot of N shows the presence of a predominant tip along the first diagonal.

VII. CONCLUSION
We have developed a model to describe a non-interacting DQD system regardless its geometry, either in series or in parallel, and derived its electrical and thermoelectrical properties. In the case of single-level dots, the expressions for the non-equilibrium Green functions and electrical current are derived exactly. In the case of multi-level 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 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 non-interacting DQD system is made in Ref. 78 following the theoretical approach developed in Refs. 79 and 80. A direct and essential extension of this work consists in taking into account the spin degrees of freedom and the Coulomb interactions in the dots.
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 the 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), (c) and (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 .

ML-DQD
For the three energy level ML-DQD considered in the Sec. V 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 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. V B to determine the equienergetic curves at equilibrium: E λ dots = µ. The results obtained numerically are displayed in Figs. 13(b)-(d)-(f) at µ = 0 for different values of the interdot coupling 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 (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. V A and V B.