# A fully analytical equation of motion approach for the double quantum dot in
the Coulomb blockade regime

Nahual Sobrino<EMAIL_ADDRESS>Nano-Bio Spectroscopy Group and
European Theoretical Spectroscopy Facility (ETSF), Departamento de Polímeros y
Materiales Avanzados: Física, Química y Tecnología, Universidad del País Vasco
UPV/EHU, Avenida de Tolosa 72, E-20018 San Sebastián, Spain  David Jacob
Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy Facility
(ETSF), Departamento de Polímeros y Materiales Avanzados: Física, Química y
Tecnología, Universidad del País Vasco UPV/EHU, Avenida de Tolosa 72, E-20018
San Sebastián, Spain  IKERBASQUE, Basque Foundation for Science, Plaza Euskadi
5, E-48009 Bilbao, Spain  Stefan Kurth  Nano-Bio Spectroscopy Group and
European Theoretical Spectroscopy Facility (ETSF), Departamento de Polímeros y
Materiales Avanzados: Física, Química y Tecnología, Universidad del País Vasco
UPV/EHU, Avenida de Tolosa 72, E-20018 San Sebastián, Spain  IKERBASQUE,
Basque Foundation for Science, Plaza Euskadi 5, E-48009 Bilbao, Spain
Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4,
E-20018 San Sebastián, Spain

###### Abstract

A fully analytical approach based on the equation of motion technique to
investigate the spectral properties and orbital occupations in an interacting
double quantum dot in equilibrium is presented. By solving a linear system for
the density correlators analytically, an explicit expression for the one body
Green’s function in terms of local occupations, intra- and inter-dot Coulomb
interactions, and the model parameters is derived. In the uncontacted limit,
the results coincide with those obtained from the grand canonical ensemble.
The analytical results compare favorably with numerical results obtained with
the Non-Crossing Approximation and the Hierarchical Equation of Motion methods
accurately reproducing peak positions and spectral weight distributions in the
Coulomb blockade regime.

##  I Introduction

In the last decades Quantum dots (QDs) have attracted an enormous amount of
attention due to their potential applications in quantum computing,
nanoelectronics, and as ideal systems to explore fundamental quantum
phenomena[1, 2, 3, 4, 5, 6, 7, 8, 9]. QDs are essentially nanostructures that
confine electrons within a small region of space, leading to the quantization
of the electronic motion in all dimensons. They can be realized in manifold
ways, e.g., by electrostatic confinement of a two-dimensional electron gas
[1], metallic or semiconductor nanoparticles [10, 11], carbon nanotubes [12,
13], or the buckminster fullerene (C60) [14].

A QD coupled to electrodes can be effectively described by an Anderson
impurity model (AIM)[15], which has been a cornerstone for understanding
strongly correlated electronic systems. The AIM offers an ideal playground to
explore electron-electron interactions, shedding light on phenomena such as
the Kondo effect and magnetic impurity behavior in metals, significantly
advancing our understanding of many-body physics[16]. These insights have
paved the way for exploring more complex systems, such as double quantum dots
(DQDs), which offer a unique platform for investigating electron-electron
interactions, quantum coherence, and transport properties in low-dimensional
systems. The coupling between QDs and the presence of interdot interactions
introduce complex and intriguing physics that have been extensively studied in
and out of equilibrium both theoretically and experimentally [17, 18, 19, 20,
21, 22, 23, 24].

Various computational approaches have been developed for solving the AIM and
its extensions, among them numerical renormalization group (NRG), quantum
Monte Carlo (QMC), non-crossing approximation (NCA) and hierarchical equations
of motion (HEOM)[25, 26, 27, 28, 29, 30, 31, 32, 33]. On the other hand, the
equation of motion (EOM) approach has emerged as a powerful tool for
investigating Green’s functions in Hubbard and Anderson impurity models[34,
35, 36, 37, 38]. The EOM approach provides a systematic way to derive the
equations governing the dynamics of the Green’s functions (GFs), allowing for
a detailed analysis of the spectral properties and electronic correlations.
This approach has been successfully applied to study the spectral and
transport properties of single and multiple quantum dot systems under various
interaction regimes[39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 46, 51].
However, the complexity of the equations due to inter-dot interactions and
coupling to electronic reservoirs, required different approximations and
numerical techniques to accurately describe correlated electronic states and
transport phenomena.

Despite the progress made using numerical techniques, an analytical approach
offers the distinct advantage of providing explicit functional dependencies.
Understanding these dependencies is crucial for designing and optimizing
quantum dot-based devices, particularly in regimes where strong electron-
electron interactions play a dominant role. Furthermore, analytical solutions
offer a computationally efficient way to study additional properties and
explore new physical regimes, guiding experimental efforts and advancing the
development of quantum technologies.

The rest of the paper is organized as follows: In Sec. II, we first introduce
the DQD Hamiltonian coupled to reservoirs, and explain in detail our analytic
solution for the DQD via the EOM approach. In Sec. III, we provide a
comprehensive comparison of the analytical results with numerical methods,
discussing the accuracy and reliability of the analytical EOM approach. The
conclusions are presented in Sec. IV. Finally, we provide details of the
derivations in three appendices: App. A details our approximate treatment of
the leads, App. B gives details for the different GFs while App. C presents
the linear system for the density correlators and its analytical solution.

##  II Model and Method

###  II.1 Double quantum dot model

We consider a parallel double quantum dot (DQD) system, where each quantum dot
(QD) $i=1,2$ is attached to an individual reservoir (thus also labeled by $i$)
of non-interacting electrons. The Hamiltonian of the system has the form

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\sum_{i}v_{i}\,\hat{n}_{i}+\sum_{i}U_{i}\,\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}+U_{12}\,\hat{n}_{1}\hat{n}_{2}$
$\displaystyle+\sum_{ik\sigma}\epsilon_{ik}\,\hat{c}^{\dagger}_{ik\sigma}\hat{c}_{ik\sigma}+\sum_{ik\sigma}\left(V_{i,k}\,\hat{c}^{\dagger}_{ik\sigma}\hat{d}_{i\sigma}+\text{H.c.}\right)$
(1)

where $\hat{d}_{i\sigma}$ ($\hat{d}_{i\sigma}^{\dagger}$) is the annihilation
(creation) operator for an electron with spin $\sigma$ on QD $i$, and
$\hat{c}_{ik\sigma}$ ($\hat{c}_{ik\sigma}^{\dagger}$) is the annihilation
(creation) operator for an electron state $k,\sigma$ in reservoir $i$.
$\hat{n}_{i}=\sum_{\sigma}\hat{n}_{i\sigma}$ and
$\hat{n}_{i\sigma}=\hat{d}^{\dagger}_{i\sigma}\hat{d}_{i\sigma}$ are the total
and the spin resolved occupation operators for QD $i$, respectively.

The first three terms in Eq. 1 describe the DQD with on-site energies (or
gates) $v_{i}$ and subject to intra- and inter-dot Coulomb repulsion $U_{i}$
and $U_{12}$, respectively. The last two terms in Eq. 1 describe the single-
particle eigenstates of the reservoirs with energies $\epsilon_{ik}$ and the
coupling to the corresponding QD with hopping $V_{i,k}$.

###  II.2 Analytical Equation of Motion approach

We now apply the EOM approach to obtain an analytic expression for the single-
particle GF solely in terms of the occupations $\langle\hat{n}_{i}\rangle$ of
each QD $i$ and the electron addition and removal energies. The derived
expression is valid for the weak-coupling regime, i.e. when the coupling to
the reservoirs is small compared to the temperature, and becomes exact in the
limit of vanishing coupling, in contrast to other approaches which usually
rely on some kind of mean-field decoupling [45, 46]. Our approach consists of
two steps: First the EOMs are solved recursively, leading to expressions for
both the single-particle as well as higher-order GFs in terms of density
correlators
$\langle\hat{n}_{i_{1}\sigma_{1}}\hat{n}_{i_{2}\sigma_{2}}\ldots\rangle$ where
the brackets $\braket{\dots}$ indicate the thermal average, i.e.,
$\braket{\hat{\mathcal{A}}}={\rm
Tr}[e^{-\beta\hat{\mathcal{H}}}\hat{\mathcal{A}}]/{\rm
Tr}[e^{-\beta\hat{\mathcal{H}}}]$ with the inverse temperature $\beta\equiv
1/T$. From the equations for the GFs, in a second step we obtain a set of
linear equations for the correlators, which in the case of the DQD can be
solved analytically. The solution allows to express the correlators and thus
the GF in terms of the dot occupations $\langle\hat{n}_{i}\rangle$.

The retarded single-particle GF of the DQD system given in Eq. 1 is defined in
the time domain as

$\displaystyle
G_{i\sigma}^{r}(t)=-i\theta(t)\braket{\\{\hat{d}_{i\sigma}(t),\hat{d}_{i\sigma}^{\dagger}(0)\\}}$
(2)

where $\theta(t)$ is the Heaviside step function , the curly brackets denote
the anticommutator, and for any operator $\cal{\hat{A}}$ the corresponding
operator in the Heisenberg picture is
$\hat{\cal{A}}(t)=\exp(i{\cal{\hat{H}}}t)\hat{\cal{A}}\exp(-i{\cal{\hat{H}}}t)$.
The Fourier transform of Eq. 2

$\displaystyle
G_{i\sigma}^{r}(\omega)\equiv\braket{\braket{\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\int_{-\infty}^{\infty}G_{i\sigma}^{r}(t)\,e^{i\omega
t}dt$ (3)

can be directly related to the orbital occupation via
($\int\equiv\int\frac{d\omega}{2\pi}$ in the following)

$\displaystyle\braket{\hat{n}_{i\sigma}}=-\int
f(\omega)\,\text{Im}\left(G_{i\sigma}^{r}(\omega)\right).$ (4)

where $f(\omega)=[1+e^{\beta\omega}]^{-1}$ is the Fermi function.

In order to proceed, it is convenient to introduce a more general GF that
represents the GF generated at all orders of the EOM approach:

$\displaystyle
G_{B}^{r}(t)=-i\theta(t)\braket{\\{\hat{B}(t),\hat{d}_{i\sigma}^{\dagger}(0)\\}}$
(5)

where the operator $\hat{B}(t)$ is a product of occupation operators
multiplied by either a single dot annihilation operator $\hat{d}_{i\sigma}$ or
a lead annihilation operator $\hat{c}_{ik\sigma}$ in the Heisenberg picture.
The EOM of its Fourier transform
$G_{B}^{r}(\omega)\equiv\braket{\braket{\hat{B}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
in the frequency domain can be written as [34, 52, 36]

$\displaystyle\omega^{+}\braket{\braket{\hat{B}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\braket{\\{\hat{B},\hat{d}^{\dagger}_{i\sigma}\\}}+\braket{\braket{[\hat{B},\hat{\mathcal{H}}]:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(6)

where $\omega^{+}\equiv\omega+i\eta$ with $\eta\rightarrow 0^{+}$ is the
energy $\omega$ shifted infinitesimally to the complex plane. We define the
general $N+1$ particle dot GF as
$G_{Dd_{i\sigma}}^{r}(\omega)\equiv\braket{\braket{\hat{D}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
where
$\hat{D}=\hat{n}_{i_{1}\sigma_{1}}\hat{n}_{i_{2}\sigma_{2}}\ldots\hat{n}_{i_{N}\sigma_{N}}$
with $(i_{k},\sigma_{k})\neq(i_{q},\sigma_{q})$ for all $k\neq q$ and
$(i_{k},\sigma_{k})\neq(i\sigma)$ for all $k\in\\{1,\ldots,N\\}$ and
$\hat{D}=\hat{\mathbbm{1}}$ is the unit operator if $N=0$.

In the following we use the EOM to derive analytical expressions for the dot
equilibrium GFs $G^{r}_{Dd_{i\sigma}}(\omega)$ of the DQD system attached to
two reservoirs solely in terms of the local occupations
$\braket{\hat{n}_{i\sigma}}$ and the removal and addition energies of the DQD.
Our starting point is the EOM for the retarded single-particle GF

$\displaystyle\omega^{+}\,\braket{\braket{\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=$
$\displaystyle\braket{\\{\hat{d}_{i\sigma},\hat{d}^{\dagger}_{i\sigma}\\}}+\braket{\braket{[\hat{d}_{i\sigma},\hat{\mathcal{H}}]:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(7)

The first term of the r.h.s. of Eq. 7 directly follows from the anti-
commutation rules for creation and annihilation operators:
$\braket{\\{\hat{d}_{i\sigma},\hat{d}^{\dagger}_{i\sigma}\\}}=1$. The
commutator in the second term of the r.h.s. of Eq. 7 is found to be

$\displaystyle[\hat{d}_{i\sigma},\hat{\mathcal{H}}]=\left(v_{i}+U_{i}\,\hat{n}_{i\bar{\sigma}}+U_{12}\,\hat{n}_{\bar{i}}\right)\hat{d}_{i\sigma}+\sum_{k}V^{*}_{i,k}\,\hat{c}_{ik\sigma}.$
(8)

Hence the second term of the r.h.s. of Eq. 7 generates a contribution to the
single-particle GF, two higher order GFs and a mixed dot-lead GF according to:

$\displaystyle(\omega^{+}-v_{i})\,G_{i\sigma}^{r}(\omega)=1+U_{i}\braket{\braket{\hat{n}_{i\bar{\sigma}}\,\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle+U_{12}\braket{\braket{\hat{n}_{\bar{i}}\,\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}+\sum_{k}V^{*}_{i,k}\braket{\braket{\hat{c}_{ik\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(9)

The EOM of the mixed dot-lead GF is obtained from Eq. 6 with
$\hat{B}=\hat{c}_{ik\sigma}$

$\displaystyle\omega^{+}\braket{\braket{\hat{c}_{ik\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\braket{\braket{[\hat{c}_{ik\sigma},\hat{\mathcal{H}}]:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(10)

where the commutator between the annihilation operator of the leads and the
Hamiltonian is
$[\hat{c}_{ik\sigma},\hat{\mathcal{H}}]=\epsilon_{ik}\hat{c}_{ik\sigma}+V_{i,k}\hat{d}_{i\sigma}$.
Then, the last term of the r.h.s of Eq. 9 becomes

$\displaystyle\sum_{k}V^{*}_{i,k}\braket{\braket{\hat{c}_{ik\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\sum_{k}\frac{V_{i,k}V_{i,k}^{*}}{\omega^{+}-\epsilon_{ik}}\,G_{i\sigma}^{r}(\omega).$
(11)

The prefactor multiplying the GF $G_{i\sigma}^{r}$ on the r.h.s. of this
equation is known as hybridization function or embedding self-energy for QD
$i$, defined by

$\Delta_{i}(\omega)\equiv\sum_{k}\frac{V_{i,k}^{\ast}\,V_{i,k}}{\omega^{+}-\epsilon_{ik}}$
(12)

The hybridization function $\Delta_{i}(\omega)$ describes the effect of the
coupling of QD $i$ to its reservoir. More specifically, its real part yields
the renormalization of the on-site energy $v_{i}$, while its imaginary part
describes the life-time broadening of the QD level due to the coupling to the
reservoir.

Taking into account the definition of the hybridization function Eq. 12, the
single-particle GF can be written as

$\displaystyle(\omega-
v_{i}-\Delta_{i}(\omega))\,G_{i\sigma}^{r}(\omega)=1+U_{i}\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle+U_{12}\braket{\braket{\hat{n}_{\bar{i}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(13)

Setting up the EOMs for the higher order GFs on the r.h.s. of (13) leads, on
the one hand, to the generation of GFs of yet higher order of the form
$\braket{\braket{\hat{n}_{i^{\prime\prime}\sigma^{\prime\prime}}\hat{n}_{i^{\prime}\sigma^{\prime}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$.
In each subsequent order the number of occupation operators in the generated
higher order GFs thus increases by one. On the other hand, it generates a
mixed dot-lead GF of the same order
$\braket{\braket{\hat{n}_{i^{\prime}\sigma^{\prime}}\hat{c}_{ik\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$.
In order to truncate the hierarchy of the the mixed GF, we now assume that the
local density operator and the Hamiltonian commute, i.e.,
$[\hat{n}_{i\sigma},\hat{\mathcal{H}}]\approx 0$. The GFs calculated within
this approximation are accurate in the CB regime, when the temperature $kT$ is
large compared to the broadening $-{\rm Im}\Delta_{i}$ by the reservoirs. This
approximation can also be obtained by neglecting certain GFs which contain one
or more lead operators (for a detailed derivation see Appendix A).

Under this approximation, the mixed higher-order GFs are related to the same-
order dot GFs via

$\displaystyle\braket{\braket{\hat{D}\hat{c}_{ik\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle\approx\braket{\braket{\hat{D}[\hat{c}_{ik\sigma},\hat{\mathcal{H}}]:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle=\frac{V_{i,k}}{\omega^{+}-\epsilon_{ik}}\braket{\braket{\hat{D}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}\;.$
(14)

Therefore, the one-body dot GF Eq. 3 can be expressed completely in terms of
higher-order dot GFs $G_{Dd_{i\sigma}}^{r}(\omega)$. Hence, for a finite
number of single-particle levels the EOM can be closed at finite order. The
EOM approach can thus be solved by recursive evaluation of the EOMs for
successive orders of GFs. We can easily see that for the DQD the EOMs can be
closed at the fourth order $G_{Dd_{i\sigma}}^{r}(\omega)$ where
$\hat{D}=\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}$,
since the additional operator $\hat{n}_{i^{\prime}\sigma^{\prime}}$ generated
by the interaction part of the EOM must be equal to either one of the
operators $\hat{n}_{i\bar{\sigma}}$, $\hat{n}_{\bar{i}\sigma}$ or
$\hat{n}_{\bar{i}\bar{\sigma}}$ in $\hat{D}$. This yields the four-particle GF

$\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}}{\omega-
v_{i}-U_{i}-2U_{12}-\Delta_{i}(\omega)}$ (15)

As shown in App. B, replacing the expressions of the higher-order GFs
(expressed in terms of the correlators Eq. 19) recursively into the lower-
order GFs finally allows to express the one-body GF (3) as

$\displaystyle G_{i\sigma}^{r}(\omega)=\sum_{j=1}^{6}\frac{r_{i,j}}{\omega-
p_{i,j}-\Delta_{i}(\omega)}$ (16)

where in order to simplify the notation, we have introduced the following
abbreviations for the addition and removal energies that define the poles of
the GF (16):

$\displaystyle p_{i,1}=v_{i}$ (17a) $\displaystyle p_{i,2}=v_{i}+U_{i}$ (17b)
$\displaystyle p_{i,3}=v_{i}+U_{i}+U_{12}$ (17c) $\displaystyle
p_{i,4}=v_{i}+U_{i}+2\,U_{12}$ (17d) $\displaystyle p_{i,5}=v_{i}+U_{12}$
(17e) $\displaystyle p_{i,6}=v_{i}+2\,U_{12}$ (17f)

Importantly, the residues $r_{i,j}$ in Eq. 16 can be expressed by linear
combinations of the one, two and three-body correlators as follows:

$\displaystyle r_{i,1}$ $\displaystyle=$ $\displaystyle
1-\braket{\hat{n}_{i\sigma}}-2\braket{\hat{n}_{\bar{i}\sigma}}+\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}+2\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}$
(18a)
$\displaystyle-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
$\displaystyle r_{i,2}$ $\displaystyle=$
$\displaystyle\braket{\hat{n}_{i\sigma}}-2\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}+\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
(18b) $\displaystyle r_{i,3}$ $\displaystyle=$ $\displaystyle
2(\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}})$
(18c) $\displaystyle r_{i,4}$ $\displaystyle=$
$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
(18d) $\displaystyle r_{i,5}$ $\displaystyle=$ $\displaystyle
2(\braket{\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}$
(18e)
$\displaystyle+\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}})\hskip
21.52771pt$ $\displaystyle r_{i,6}$ $\displaystyle=$
$\displaystyle\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
(18f)

Eq. 16 provides the structure of the one body GF as a sum of single poles,
broadened by the imaginary part of the hybridization function $\Delta_{i}$. In
the wide-band limit (WBL) the hybridization function becomes constant and
purely imaginary, $\Delta_{i}=-i\gamma/2$. We emphasize that also the higher-
order dot GFs have a similar structure as Eq. 16 but in this case only a
subset of the poles of Eq. 17 contribute, see Appendix B. The imaginary
contribution in the denominators is obtained under the approximation
$[\hat{n}_{i\sigma},\hat{\mathcal{H}}]=0$ and it formally justifies the usual
phenomenological procedure in which $\gamma$ is introduced as a decay
parameter in order to broaden the poles of the spectral function for different
impurity systems [53, 54, 55].

The second step in our analytical EOM approach rests on the observation that
one can obtain the static $N$-body correlator from the GF via [34, 56]

$\displaystyle\braket{\hat{D}\,\hat{n}_{i\sigma}}=-\int
f(\omega)\,\text{Im}\left(G_{Dd_{i\sigma}}^{r}(\omega)\right).$ (19)

Since the GFs $G_{Dd_{i\sigma}}^{r}(\omega)$ only depend linearly on the
correlators, application of Eq. 19 for all GFs then leads to a linear system
for the correlators. Importantly, the fourth-order correlator
$\braket{\hat{n}_{1\sigma}\hat{n}_{1\bar{\sigma}}\hat{n}_{2\sigma}\hat{n}_{2\bar{\sigma}}}$
does not explicitly appear in the resulting linear system which thus in
principle consists of 15 equations for 15 independent correlators.

The size of the linear system can be reduced due to the spin symmetry in our
model (the onsite energies $v_{i}$ are independent of spin) which reduces the
number of independent correlators from 15 to 7 (see Appendix C). Incidentally,
in the method of rate equations [57, 58, 59] one solves for the same number of
independent variables which in this case are the probabilities for a certain
distribution of occupation numbers in the double dot. In fact, we expect the
rate equation approach to be physically equivalent to the EOM approach at our
level of approximate treatment of the coupling to the leads.

Another important simplification arises from the fact that all GFs have the
analytic structure of a sum of single poles in the complex plane and therefore
each integral arising on the r.h.s. of Eq. 19 has the same structure and can
be evaluated analytically [60, 47]

$\displaystyle\phi(p)$ $\displaystyle=\int
f(\omega)\frac{\gamma}{(\omega-p)^{2}+\frac{\gamma^{2}}{4}}$
$\displaystyle=\frac{1}{2}-\frac{1}{\pi}\psi\left(\frac{1}{2}+\frac{\gamma/2+ip}{2\pi
T}\right)$ (20)

where $\psi(z)=\frac{d\log(\Gamma(z))}{dz}$ is the digamma function with
general complex argument $z$, and $\Gamma(z)$ is the gamma function. For
instance, the two-body correlator function
$\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}$, can be
expressed as

$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}=-\int
f(\omega)\text{Im}\left(\braket{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\bar{\sigma}}:\hat{d}^{\dagger}_{i\bar{\sigma}}}}\right)$
$\displaystyle=\phi(p_{i,5})\left(\braket{\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}+\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\right)$
$\displaystyle\;\;+\phi(p_{i,6})\left(\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\right)$
$\displaystyle\;\;+\phi(p_{i,3})\left(\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\right)$
$\displaystyle\;\;+\phi(p_{i,4})\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\;.$
(21)

In other words, all coefficients in the linear system for the correlators can
be computed analytically. Of course, once these coefficients are known, the
resulting linear system can readily be solved numerically. However, in the
case of our double dot one can even solve the whole linear system
analytically. To this end, for the quantum dot with label $i$ we first solve
the subsystem for the higher-order correlators
$\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}$,
$\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}$,
$\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$,
and $\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
and express them in terms of the occupations
$\braket{\hat{n}_{i\sigma}}=n_{i}/2$, and
$\braket{\hat{n}_{\bar{i}\sigma}}=n_{\bar{i}}/2$. Subsequently, we solve the
resulting system for the two local occupations (see Appendix C for the full
derivation).

The higher-order correlators are expressed in terms of the occupations as

$\displaystyle\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}=\tau_{i,1}\braket{\hat{n}_{\bar{i}\sigma}}$
(22a)
$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}=\tau_{i,2}\braket{\hat{n}_{\bar{i}\sigma}}$
(22b)
$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}=\tau_{i,3}\braket{\hat{n}_{\bar{i}\sigma}}$
(22c)
$\displaystyle\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}=\tau_{i,4}\braket{\hat{n}_{\bar{i}\sigma}}$
(22d)

where the parameters $\tau_{i,j}$ are defined in Eq. 57 and are fully given in
terms of the model parameters. Substituting Eq. 22 into Eq. 18 allows to
express the residues of the single-particle GF in terms of the orbital
occupations

$\displaystyle r_{i,1}=$ $\displaystyle
1-n_{i}/2+\left(\tau_{i,1}+2\tau_{i,2}-\tau_{i,3}-2\right)n_{\bar{i}}/2$ (23a)
$\displaystyle r_{i,2}=$ $\displaystyle
n_{i}/2+\left(\tau_{i,3}-2\tau_{i,2}\right)n_{\bar{i}}/2$ (23b) $\displaystyle
r_{i,3}=$ $\displaystyle\left(\tau_{i,2}-\tau_{i,3}\right)n_{\bar{i}}$ (23c)
$\displaystyle r_{i,4}=$ $\displaystyle\tau_{i,3}n_{\bar{i}}/2$ (23d)
$\displaystyle r_{i,5}=$
$\displaystyle\left(1-\tau_{i,1}-\tau_{i,2}+\tau_{i,3}\right)n_{\bar{i}}$
(23e) $\displaystyle r_{i,6}=$
$\displaystyle\left(\tau_{i,1}-\tau_{i,3}\right)n_{\bar{i}}/2$ (23f)

Combining Eqs. 4, 16 and 23 one then arrives at a $2\times 2$ linear system
for the local occupations which is readily solved analytically, see Appendix
C. The resulting local occupations are then given by:

$\displaystyle
n_{i}=2\frac{\phi(\varepsilon_{i})\eta_{\bar{i}\bar{i}}-\phi(\varepsilon_{\bar{i}})\eta_{i\bar{i}}}{\eta_{11}\eta_{22}-\eta_{12}\eta_{21}}$
(24)

where

$\displaystyle\eta_{ii}=$ $\displaystyle 1+\phi(p_{i,1})-\phi(p_{i,2})$ (25a)
$\displaystyle\eta_{i\bar{i}}=$
$\displaystyle\phi(p_{i,1})\left(2+\tau_{i,3}-\tau_{i,1}-2\tau_{i,2}\right)+\phi(p_{i,2})\left(2\tau_{i,2}-\tau_{i,4}\right)$
$\displaystyle+$ $\displaystyle
2\phi(p_{i,3})\left(\tau_{i,3}-\tau_{i,2}\right)-\tau_{i,3}\phi(p_{i,4})$
$\displaystyle+$ $\displaystyle
2\phi(p_{i,5})\left(\tau_{i,1}+\tau_{i,2}-\tau_{i,3}-1\right)+\phi(p_{i,6})\left(\tau_{i,3}-\tau_{i,1}\right)$
(25b)

Eqs. 16, 23 and 24 represent the main results of this work, providing
analytical expressions for the evaluation of the Green’s functions and the
orbital occupations as functions of the model parameters.

In the uncontacted limit $\gamma\to 0$, our method involves no approximations
and treats interactions exactly. The imaginary part of the resulting GF (Eq.
16) is proportional to a sum of delta functions at frequencies coinciding with
the poles $p_{i,j}$ which allows to express the correlators in terms of Fermi
functions, $\phi(p)=f(p)$. The results in this limit are in complete agreement
with the grand canonical ensemble (GCE) as it should be.

A common alternative approach to treat the coupling dependence consists in
truncating the mixed GFs hierarchy at the level of the two-particle GF,
considering the Hubbard-I [35] decoupling
$\braket{\braket{\hat{n}_{i^{\prime}\sigma^{\prime}}\hat{c}_{ik\sigma}:d^{\dagger}_{i\sigma}}}\approx\braket{\hat{n}_{i^{\prime}\sigma^{\prime}}}\braket{\braket{\hat{c}_{ik\sigma}:d^{\dagger}_{i\sigma}}}$
and treating the inter-dot Coulomb interaction within the Hartree
approximation. This method (denoted as EOM-H in the following) has already
been explored in the literature [46] and the resulting one-particle GF is

$\displaystyle G_{i\sigma}^{r}(\omega)$
$\displaystyle=\frac{1-n_{i}/2}{\omega-
p_{i,1}-U_{12}n_{\bar{i}}-\Delta_{i}\left(1+\frac{U_{i}\,n_{i}/2}{\omega-
p_{i,2}}\right)}$ $\displaystyle+\frac{n_{i}/2}{\omega-
p_{i,2}-U_{12}n_{\bar{i}}-\Delta_{i}\left(1-\frac{U_{i}(1-n_{i}/2)}{\omega-
p_{i,1}}\right)}$ (26)

Note that since the spectral function related to Eq. 26 can not be written as
a sum of Lorentzian functions, the resulting integrals for obtaining the
occupations have to be solved numerically. Furthermore, the resulting system
of equations is nonlinear in the occupations which cannot be solved
analytically and requires a numerical approach.

##  III Results

In this section, we apply our analytical EOM approach to compute spectral
functions and orbital occupations for the DQD in different regimes defined by
the interaction parameters (c.f. Ref. 61). We benchmark the results with two
well-known numerical methods, namely the non-crossing approximation (NCA) [29,
30] and the hierachical equations of motion (HEOM) approach [62, 63, 64]. A
Python code implementing our EOM approach for the DQD is publicly available at
GitHub [65].

Figure 1:  Green’s function residues $r_{1,j}$ of Eq. 23 as a function of the
gate potentials for different combinations of the interaction parameters at
$T=0.05$ (upper panels) and $T=0.5$ (lower panels) in the limit $\gamma\to 0$.
The Coulomb interactions are: left) $U_{1}=2$ and $U_{12}=0.5$, center)
$U_{1}=2$ and $U_{12}=1.3$, and right) $U_{1}=2$ and $U_{12}=1.8$. Energies
given in units of $U_{2}$.

In the following we consider the three different regimes of interaction
parameters $U_{1}$, $U_{2}$ and $U_{12}$ identified in Ref. 61: In Regime I
the inter-dot interaction is smaller than either of the intra-dot
interactions. Without loss of generality we may also assume
$U_{2}\leq{U_{1}}$, and hence $U_{12}\leq{U_{2}}\leq{U_{1}}$. Regime I is the
most common (or natural) case, since in physical systems usually the inter-dot
(or inter-orbital) interaction will be smaller than the intra-dot (or intra-
orbital) ones. On the other hand, Regimes II ($(U_{1}+U_{2})/2>U_{12}>U_{2}$)
and III ($U_{12}>(U_{1}+U_{2})/2>U_{2}$) where $U_{12}$ is larger than at
least one of the on-site interactions $U_{i}$ could be relevant in specific
setups or for effective model Hamiltonians.

In order to better understand the structure of the GF given by Eq. 16, in Fig.
1 we plot the residues given by Eq. 23 for QD $i=1$ as functions of the gate
potentials $v_{i}$ of both QDs for each of the three aforementioned regimes,
together with so-called stability diagrams structure(shown as dashed grey
lines) in the limit of low coupling and both low (upper panels) and higher
(lower panels) temperatures (on the scale of the interactions). A stability
diagram is a map of the dot occupation numbers $\\{n_{1},n_{2}\\}$ as a
function of the gate potentials $\\{v_{1},v_{2}\\}$. It thus shows the regions
where the occupation numbers are stable, i.e., do not fluctuate in the limit
of low temperature and weak coupling to the reservoirs. Specifically in this
limit and for the case of only density-density interactions as in our case,
the occupations in the different regions are stable at integer values. As
explained in detail in Ref. 61, each regime is characterized by a distinctive
stability diagram.

Regime I, shown in Fig. 1(a), is characterized by relatively independent
behavior of both dots and a wide central region of occupations $(1,1)$. In
Regime II, shown in Fig. 1(b), the width of the $(1,1)$ region in the
stability diagram is notably reduced, leading to increased regions of the
$(1,2)$ and $(2,1)$ occupation. Finally, in Regime III, the $(1,1)$ region of
occupations is not present anymore, and the regions $(0,1)$ and $(2,0)$ as
well as the regions $(1,0)$ and $(0,2)$ become adjacent.

We now focus on a given regime, i.e., fixed interactions $U_{1}$, $U_{2}$ and
$U_{12}$, and a given region in the $v_{1}-v_{2}$ plane. By inspection of the
upper panels of Fig. 1 ( low coupling and low temperature), we realize that
for each given region there are two situations possible: (i) there exists only
one non-vanishing residue with value 1 or (ii) there exist two non-vanishing
residues both with value 0.5, i.e., the local spectral function (on dot 1) has
either one or two poles. In the first case the single pole is at an electron
removal or electron addition energy for an electron at the given site. In the
case of two poles, the corresponding site is half filled and one can both add
and remove an electron to/from that site. Note that both poles have equal
spectral weight, i.e., equal residues, because the dot Hamiltonian only
contains density-density interactions. As the temperature is increased,
eigenstates different from the ground state start to contribute to the
ensemble of the dot, and thus more than two poles may contribute to the local
spectral function, as can be seen in the lower panels of Fig. 1 (higher
temperature).

Figure 2: Spectral function of the double quantum dot from Eqs. 16 and 23 for
$T=1/6$, $\gamma=1/30$ at the particle-hole symmetric point
$v_{i}=-U_{i}/2-U_{12}$. (a) $U_{1}=2/3$ and (b) $U_{1}=0$. Energies in units
of $U_{2}$.

In Fig. 2, the total spectral function
$A(\omega)=-\text{Im}[G_{1}^{r}(\omega)+G_{2}^{r}(\omega)]/\pi$ obtained from
the analytical EOM, is plotted as a function of the inter-Coulomb repulsion
$U_{12}$ at the particle-hole symmetric point $v_{i}=-U_{i}/2-U_{12}$, with
$T/\gamma=5$, $U_{2}=1$, and (a) $U_{1}=2/3$ and (b) $U_{1}=0$. The main
structure of the spectral function consists of a series of intersecting
vertical and diagonal lines. The vertical lines are located at $\omega=\pm
U_{i}/2$ while the diagonals follow from $\pm\omega=\pm U_{i}/2-U_{12}$
consistent with the position of the poles Eq. 17 for $v_{i}=-U_{i}/2-U_{12}$.

An important observation is that for small $U_{12}$ values, almost all the
spectral weight is carried by the vertical lines. As $U_{12}$ increases, the
spectral weight shifts to the diagonal lines. This occurs at the transitions
to regime II (at (panel a) $U_{12}=2/3$ and to regime III (at (panel a)
$U_{12}=5/6$ and (panel b) $U_{12}=1/2$). The lines with less spectral weights
are due to excited states and their weight decreases in the limit of both
$T\to 0$ and $\gamma\to 0$ (not shown).

Figure 3:  Spectral function of the double quantum dot at $T=1/6$,
$\gamma=1/30$ and at the particle-hole symmetric point
$v_{i}=-U_{i}/2-U_{12}$. In (a) $U_{1}=2/3$, $U_{12}=7/6$, (b) $U_{1}=2/3$,
$U_{12}=1$, (c) $U_{1}=0$, $U_{12}=7/6$, and (d) $U_{1}=0$, $U_{12}=1/6$. The
different lines correspond to the EOM Hartree approximation (EOM-H) (divided
by 2 for representation clarity), the analytical EOM of Eqs. 16 and 23, and
the numerical results from NCA and HEOM. Energies in units of $U_{2}$.

In Fig. 3, we compare our results for the total spectral function (EOM) for
those values of $U_{12}$ marked by the green dashed lines in Fig. 2 against
three other methods: the Hubbard I + Hartree decoupling scheme of Eq. 26
(EOM-H), NCA, and HEOM for tier level $L=3$. At small $U_{12}$ values, the
EOM-H approach captures the general distribution of the spectral functions
(corresponding to the vertical lines in Fig. 2) in comparison to NCA and HEOM.
However, as $U_{12}$ increases, EOM-H fails to capture the correct spectral
weights of the different peaks, finding that the vertical lines still carry
all the spectral weight (instead of the diagonal ones). As shown, the
analytical EOM accurately reproduces the position of all the peaks and
correctly captures the spectral weight distributions for all values of
$U_{12}$.

Figure 4: Local occupations (left) and spectral functions (right) as function
of the average gate level $v=(v_{1}+v_{2})/2$ for $U_{12}=0.5$, $T=0.1$, and
$\gamma=0.05$. The other interaction and parameters are (a) $U_{1}=0.75$ and
$\delta v=0$, (b) $U_{1}=0.75$ and $\delta v=0.5$, (c) $U_{1}=0.5$ and $\delta
v=0$, and (d) $U_{1}=0.5$ and $\delta v=0.5$. Energies in units of $U_{2}$.

Finally, in the left columns of Fig. 4 we present a comparison of the local
occupations as a function of the average gate level $v=(v_{1}+v_{2})/2$
between our analytical EOM approach, the EOM-H approximation and the HEOM
method for different interactions and $\delta v=v_{1}-v_{2}$, at a fixed
$T=2\gamma=0.1$. In panel a) with $U_{1}=0.75$, the local occupations are
degenerate up to half occupation at $v\sim-0.5$. Then, $n_{1}$ evolves faster
to full occupation due to the smaller intra-dot Coulomb repulsion. In panel b)
where a finite energy difference is applied between the gate levels, $\delta
v=v_{1}-v_{2}=0.5$, the occupation of site one is shifted to lower gates, with
two plateaus close to half filling being the main features of the densities.
In panel c) where $U_{1}=0.5$ and $\delta v=0$, both $n_{i}$ evolve together
as the gate is decreased up to quarter occupation at $v\sim 0$. As in a), then
$n_{1}$ evolves faster to full occupation. Finally, in in panel d), the local
occupation of site one is shifted to lower gate levels $v$. In all the
situations studied with $T/\gamma>1$, the agreement between the EOM approach
and the HEOM is excellent, with our analytical EOM approach accurately
capturing the local step structures due to the interactions and the
corresponding slopes governed by the combination of temperature and coupling
strength. On the other hand, although the EOM-H approximation captures the
general trend of the local occupations, it fails to accurately predict some
local transitions. In particular, the EOM-H approximation does not correctly
capture the new plateaus at non-integer values of the local occupations
observed in the EOM and HEOM results. This discrepancy highlights the
limitation of the EOM-H method in dealing with the the inter-Coulomb
repulsion, not recovering correctly the GCE limit (not shown).

The right panels in Fig. 4 depict the spectral functions corresponding to each
case on the left as a function of the average gate $v$. As expected, all non-
vanishing spectral features depend linearly on $v$ due to the positions of the
poles given by $\omega=p_{i,j}$. The weight of each contribution (given by the
residues) varies as the interaction changes, but at full occupation, all the
spectral weight is carried by $r_{i,4}$ since in this case one can only remove
an electron with a removal energy given by $p_{i,4}$. In the empty dot case,
all the weight is carried by the residue $r_{i,1}$, since now one can only add
an electron with an addition energy given by $p_{i,1}$.

##  IV Conclusions

In this paper, we have introduced a fully analytical EOM approach to analyze
the spectral properties and orbital occupations of an interacting DQD system
connected to leads. A critical aspect of our approach is the systematic
closure of the hierarchy of higher-order GFs, which arises from a consistent
approximate treatment of the coupling to the leads for the GFs at all orders.
This closure is essential to map the EOMs for the GFs onto a linear system for
the density correlators, a mapping which can be applied to any number of
levels. For our double dot system it allows to derive tractable analytical
expressions for the one-body GF in terms of the local occupations and all the
model parameters. In the low coupling limit, our method correctly recovers the
results of the grand canonical ensemble, further validating the approach.

We compared our analytical results with those obtained from the NCA and the
HEOM methods. Our approach accurately reproduces the spectral peak positions
and weight distributions, demonstrating its effectiveness in capturing the
essential physics of the strongly correlated electron system in the Coulomb
blockade regime. The agreement with these established numerical techniques
underscores the reliability of our approach.

The analytical EOM approach offers substantial advantages, including explicit
functional dependencies that provide insights into the physical mechanisms
governing the system. Furthermore, the computational efficiency of our
analytical method enables a cheap exploration of different parameter regimes
of the system and the study of properties such as the transport coefficients.
Additionally, this computational efficiency may be particularly beneficial
for, e.g., density functional theory and its extensions, where knowledge of
the many-body spectral function in analytic form can be of crucial importance
for the construction of the exchange-correlation potentials by reverse
engineering.[66, 67, 68, 69, 70].

##  Appendix A Approximate treatment of coupling to leads

The purpose of the present Appendix is to treat the coupling to the leads
directly from the EOM and highlight the approximations which lead to the
contribution of the standard embedding self energy of Eq. 12 in both the
single-particle and higher order GFs (which in the WBL leads to simple
broadening of all poles of the GFs by the wide-band parameter $\gamma$ used in
Sec. II).

In Sec. II we have defined the Green function
$\langle\langle\hat{D}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
with the operator
$\hat{D}=\hat{n}_{i_{1}\sigma_{1}}\hat{n}_{i_{2}\sigma_{2}}\ldots\hat{n}_{i_{N}\sigma_{N}}$
being a product of $N$ density operators with $(i_{k}\sigma_{k})\neq(i\sigma)$
for all $k\in\\{1,\ldots,N\\}$. In addition, for $N=0$ we define
$\hat{D}=\hat{\mathbbm{1}}$ as the unit operator. We note that
$\hat{d}_{i\sigma}$ commutes with $\hat{D}$ for all $N$, i.e.,
$[\hat{D},\hat{d}_{i\sigma}]=[\hat{D},\hat{d}_{i\sigma}^{\dagger}]=0$.

After straightforward, if somewhat tedious calculations one obtains the EOM
for
$\langle\langle\hat{D}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
which reads

$\displaystyle(\omega-
v_{i})\langle\langle\hat{D}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
(27) $\displaystyle=$
$\displaystyle\langle\hat{D}\rangle+U_{i}\langle\langle\hat{D}\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{j\sigma^{\prime}}U_{ij}\langle\langle\hat{D}\hat{n}_{j\sigma^{\prime}}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{k}V_{ik\sigma}^{*}\langle\langle\hat{D}_{N}\hat{c}_{ik\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{jk\sigma^{\prime}}V_{jk\sigma^{\prime}}^{*}\langle\langle[\hat{d}_{j\sigma^{\prime}}^{\dagger},\hat{D}]\hat{d}_{i\sigma}\hat{c}_{jk\sigma^{\prime}}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{jk\sigma^{\prime}}V_{jk\sigma^{\prime}}\langle\langle\hat{c}_{jk\sigma^{\prime}}^{\dagger}[\hat{d}_{j\sigma^{\prime}},\hat{D}]\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$

Similarly, the EOM for the mixed dot-lead GF
$\langle\langle\hat{D}\hat{c}_{ik\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
reads

$\displaystyle(\omega-\varepsilon_{ik})\langle\langle\hat{D}\hat{c}_{ik\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
(28) $\displaystyle=$ $\displaystyle
V_{ik\sigma}\langle\langle\hat{D}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{jq\sigma^{\prime}}V_{jq\sigma^{\prime}}^{*}\langle\langle[\hat{D},\hat{d}_{j\sigma^{\prime}}^{\dagger}]\hat{c}_{jq\sigma^{\prime}}\hat{c}_{ik\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{jq\sigma^{\prime}}V_{jq\sigma^{\prime}}\langle\langle\hat{c}_{jq\sigma^{\prime}}^{\dagger}\hat{c}_{ik\sigma}[\hat{d}_{j\sigma^{\prime}},\hat{D}]:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$

If in Eqs. (27) and (28) we make the approximations

$\displaystyle\langle\langle[\hat{d}_{j\sigma^{\prime}}^{\dagger},\hat{D}]\hat{d}_{i\sigma}\hat{c}_{jk\sigma^{\prime}}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}\approx
0$
$\displaystyle\langle\langle\hat{c}_{jk\sigma^{\prime}}^{\dagger}[\hat{d}_{j\sigma^{\prime}},\hat{D}]\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}\approx
0$
$\displaystyle\langle\langle[\hat{D},\hat{d}_{j\sigma^{\prime}}^{\dagger}]\hat{c}_{jq\sigma^{\prime}}\hat{c}_{ik\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}\approx
0$
$\displaystyle\langle\langle\hat{c}_{jq\sigma^{\prime}}^{\dagger}\hat{c}_{ik\sigma}[\hat{d}_{j\sigma^{\prime}},\hat{D}]:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}\approx
0$ (29)

and inserting the resulting Eq. (28) into (27) we obtain

$\displaystyle(\omega-
v_{i}-\Delta_{i}(\omega))\langle\langle\hat{D}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
(30) $\displaystyle=$
$\displaystyle\langle\hat{D}\rangle+U_{i}\langle\langle\hat{D}\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$
$\displaystyle+\sum_{j\sigma^{\prime}}U_{ij}\langle\langle\hat{D}\hat{n}_{j\sigma^{\prime}}\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$

where we have used the embedding self energy $\Delta_{i}$ defined in Eq. (12)
which in the wide band limit just becomes $\Delta_{i}=-i\frac{\gamma}{2}$. We
also note that for the single-particle GF, Eq. (29) is exact and we obtain Eq.
13.

##  Appendix B Green functions

Since Eq. 1 is time-independent, the equation of motion of a general Green’s
function $\braket{\braket{\hat{B}(t):\hat{d}_{i\sigma}^{\dagger}(0)}}$ can be
expressed in the frequency domain trough the relation [52]

$\omega\braket{\braket{\hat{B}:\hat{d}_{i\sigma}^{\dagger}}}_{\omega}=\braket{\\{\hat{B},\hat{d}_{i\sigma}^{\dagger}\\}}+\braket{\braket{[\hat{B},\hat{\mathcal{H}}]:\hat{d}_{i\sigma}^{\dagger}}}_{\omega}.$
(31)

The one body GF
$\braket{\braket{\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}}}_{\omega}$ is
then derived from

$\omega\braket{\braket{\hat{d}_{i\sigma}:\hat{d}_{i\sigma}^{\dagger}}}_{\omega}=\braket{\\{\hat{d}_{i\sigma},\hat{d}_{i\sigma}^{\dagger}\\}}+\braket{\braket{[\hat{d}_{i\sigma},\hat{\mathcal{H}}]:\hat{d}_{i\sigma}^{\dagger}}}_{\omega}$
(32)

with the anti commutation relations for the creation and annihilation
operators $\\{\hat{d}_{i\sigma},\hat{d}_{j\sigma}^{\dagger}\\}=\delta_{ij}$,
$\\{\hat{d}_{i\sigma},\hat{d}_{j\sigma}\\}=\\{\hat{d}_{i\sigma}^{\dagger},d_{j\sigma}^{\dagger}\\}=0$.
The commutator of the annihilation operators in the dots of the DQD
hamiltonian follow

$\displaystyle[\hat{d}_{i\sigma},\hat{\mathcal{H}}]=\left(\varepsilon_{i}+U_{i}\hat{n}_{i\bar{\sigma}}+U_{12}\hat{n}_{\bar{i}}\right)\hat{d}_{i\sigma}+\sum_{k}V^{*}_{i,k}\hat{c}_{ik\sigma}.$
(33)

Using Eq. (30) for $\hat{D}=\hat{\mathbbm{1}}$, the equation of motion for the
single-particle GF reads

$\displaystyle\braket{\braket{\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle=\frac{1}{\omega-
v_{i}-\Delta_{i}}+\frac{U_{i}\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-\Delta_{i}}$ $\displaystyle+\frac{U_{12}}{\omega-
v_{i}-\Delta_{i}}\sum_{\sigma^{\prime}}\braket{\braket{\hat{n}_{\bar{i}\sigma^{\prime}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(34)

The two-body GF
$\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
follows from the evaluation of Eq. (30) with $\hat{D}=\hat{n}_{i\bar{\sigma}}$
which gives

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-\Delta_{i}}+$
$\displaystyle\sum_{\sigma^{\prime}}\frac{U_{12}}{\omega-
v_{i}-U_{i}-\Delta_{i}}\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma^{\prime}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
(35)

while the other two two-particle GFs are

$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{\bar{i}\sigma}}}{\omega-
v_{i}-U_{12}-\Delta_{i}}+$
$\displaystyle\frac{U_{i}\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{12}-\Delta_{i}}+\frac{U_{12}\braket{\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{12}-\Delta_{i}}$ (36)

and

$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{\bar{i}\bar{\sigma}}}}{(\omega-
v_{i}-U_{12}-\Delta_{i})}+$
$\displaystyle\frac{U_{i}\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{12}-\Delta_{i}}+\frac{U_{12}\braket{\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{12}-\Delta_{i}}$ (37)

Similarly, for the three three-body GFs we obtain

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{U_{12}\braket{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;,$ (38)
$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{U_{12}\braket{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;,$ (39)

and

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{U_{12}\braket{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;,$ (40)

Finally, the four-body GF reads

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}}{\omega-
v_{i}-U_{i}-2U_{12}-\Delta_{i}}.$ (41)

Eq. 41 is the highest order GF in the EOM hierarchy. We see that all the GFs
depend not only on the parameters of the problem but also on the various
density correlators. In order to obtain expressions for the GFs, we start by
inserting the four-particle GF of Eq. 41 into Eqs. (38)-(40) for the different
three-particle GFs. Using partial fraction decompositions along the way we
obtain

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;,$ (42)
$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;,$ (43)

and

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}=\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}\;.$ (44)

Inserting Eqs. (42)-(44) into Eqs. (35)-(37) leads to the two-particle GFs

$\displaystyle\braket{\braket{\hat{n}_{i\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle=\frac{\braket{\hat{n}_{i\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}+\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}+\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}-2\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-2U_{12}-\Delta_{i}}$ (45)
$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\sigma}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle=\frac{\braket{\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}+\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}+\frac{\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-2U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-2U_{12}-\Delta_{i}}$ (46)

and

$\displaystyle\braket{\braket{\hat{n}_{\bar{i}\bar{\sigma}}\hat{d}_{i\sigma}:\hat{d}^{\dagger}_{i\sigma}}}_{\omega}$
$\displaystyle=\frac{\braket{\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}+\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-U_{12}-\Delta_{i}}+\frac{\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}-\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-2U_{12}-\Delta_{i}}$
$\displaystyle+\frac{\braket{\hat{n}_{i\bar{\sigma}}n_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}}{\omega-
v_{i}-U_{i}-2U_{12}-\Delta_{i}}\;.$ (47)

Finally, in order to obtain the one-particle GF, we insert Eqs. (45)-(47) into
Eq. (34) to arrive at Eq. (16). In this very last step we used some of the
relations due to the spin symmetry of our problem derived in Appendix C.

##  Appendix C Linear system for density correlators and analytic solution

In the previous Appendix we have expressed the GFs in terms of density
correlators. The correlators themselves can be obtained from the GFs using Eq.
19. Thus this latter relation allows to derive a linear system for the
correlators which can easily be solved numerically. Before deriving this
linear system, however, we show that as a consequence of the spin symmetry in
our problem, i.e., the fact that the on-site energies $v_{i}$ are independent
of $\sigma$, some of the density correlators coincide. This allows to reduce
the size of the linear problem and for the case of our doube dot even to find
a fully analytic solution.

We start by using Eq. 19 on Eq. (41) with the embedding self energy in the
wide-band limit $\Delta_{i}=-i\frac{\gamma}{2}$. This leads to

$\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}=\phi(v_{i}+U_{i}+2U_{12})\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
(48)

where we used the function $\phi(p)$ defined in Eq. 20. For given dot index
$i$ the r.h.s. of Eq. 48 is independent of $\sigma$. Since the prefactor
$\phi(v_{i}+U_{i}+2U_{12})$ on the l.h.s. of the same equation is not only
strictly positive at finite temperature but also independent of $\sigma$, we
conclude that the three-body correlator
$\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}$
is independent of $\sigma$ as well.

Using this result and a similar reasoning for Eqs. (42) and (44) we find that
all two-body correlators of the form
$\braket{\hat{n}_{1\sigma}\hat{n}_{2\sigma^{\prime}}}$ are independent of both
$\sigma$ and $\sigma^{\prime}$. The spin symmetry of the two-body and the
three-body correlators can in an analogous way be used for Eq. (45) to find
that the spin densities $\braket{\hat{n}_{i\sigma}}$ on dot $i$ are
independent of $\sigma$ as expected.

If we use Eq. 19 on Eqs. (41) - (47) as well as Eq. (16) for all possible
values for the indices but do not make use of the spin symmetry we arrive at a
$14\times 15$ linear system for the one-, two-, and three-body correlators.
Instead when using the spin symmetry relations derived above, the
corresponding linear system reduces to $7\times 7$ thus leading to a
significantly reduced complexity. To be explicit, for the vector of density
correlators

${\mathbf{y}}=\left(\begin{array}[]{c}\braket{\hat{n}_{1}}\\\
\braket{\hat{n}_{2}}\\\ \braket{\hat{n}_{1\uparrow}\hat{n}_{1\downarrow}}\\\
\braket{\hat{n}_{2\uparrow}\hat{n}_{2\downarrow}}\\\
\braket{\hat{n}_{1\uparrow}\hat{n}_{2\downarrow}}\\\
\braket{\hat{n}_{1\uparrow}\hat{n}_{1\downarrow}\hat{n}_{2\uparrow}}\\\
\braket{\hat{n}_{2\uparrow}\hat{n}_{2\downarrow}\hat{n}_{1\uparrow}}\end{array}\right)$
(49)

the linear system we need to solve reads

${\mathbf{Cy}}={\mathbf{b}}$ (50)

where both the matrix $\mathbf{C}$ and the vector $\mathbf{b}$ are explicitly
given in terms of the model parameters as

$\displaystyle{\mathbf{C}}=\begin{pmatrix}1+\ell_{11}-\ell_{12}&2(\ell_{11}-\ell_{15})&0&2\ell_{15}-\ell_{11}-\ell_{16}&2\zeta_{1}&0&\nu_{1}\\\
2(\ell_{21}-\ell_{25})&1+\ell_{21}-\ell_{22}&2\ell_{25}-\ell_{21}-\ell_{26}&0&2\zeta_{2}&\nu_{2}&0\\\
-\ell_{12}&0&1&0&2(\ell_{12}-\ell_{13})&0&\ell_{13}-\ell_{12}-\ell_{14}\\\
0&-\ell_{22}&0&1&2(\ell_{22}-\ell_{23})&\ell_{23}-\ell_{22}-\ell_{24}&0\\\
0&-\ell_{15}&0&\ell_{15}-\ell_{16}&1+\ell_{15}-\ell_{13}&0&\ell_{13}+\ell_{16}-\ell_{14}-\ell_{15}\\\
0&0&0&0&-\ell_{13}&1&\ell_{13}-\ell_{14}\\\
0&0&0&0&-\ell_{23}&\ell_{23}-\ell_{24}&1\end{pmatrix}$ (51)

where we have defined

$\displaystyle\ell_{ij}=$ $\displaystyle\phi(p_{i,j})$ (52a)
$\displaystyle\zeta_{i}=$
$\displaystyle\ell_{i2}+\ell_{i5}-\ell_{i1}-\ell_{i3}$ (52b)
$\displaystyle\nu_{i}=$
$\displaystyle\ell_{i1}+2\ell_{i3}+\ell_{i6}-\ell_{i2}-\ell_{i4}-2\ell_{i5}$
(52c)

and

${\mathbf{b}}=\begin{pmatrix}\phi(p_{1,1})\\\ \phi(p_{2,1})\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\end{pmatrix}.$ (53)

While the linear system (50) can easily be solved numerically, one can also
proceed to a fully analytical solution. To this end, for a given dot index $i$
we first express the density correlators of order higher than one in terms of
the densities. This gives the smaller linear system (note that for index $i$
the correlator $\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}}$ does not
appear)

$\displaystyle\textbf{B}\begin{pmatrix}\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\\\
\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}\\\
\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}\\\
\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}\end{pmatrix}=\begin{pmatrix}0\\\
0\\\ \phi(p_{\bar{i},2})\braket{n_{\bar{i}\sigma}}\\\
\phi(p_{i,5})\braket{n_{\bar{i}\sigma}}\end{pmatrix}.$ (54)

where the matrix ${\mathbf{B}}$ is given by

$\displaystyle\textbf{B}=\begin{pmatrix}1&0&\ell_{i6}(\ell_{i4}-\ell_{i6}-1)^{-1}&0\\\
\ell_{i3}-\ell_{i4}&1&0&-\ell_{i3}\\\
0&0&1-\frac{\ell_{i6}(\ell_{\bar{i}2}+\ell_{\bar{i}4}-2\ell_{\bar{i}3})(\ell_{i4}-\ell_{i3})}{1-\ell_{i4}+\ell_{i6}}&2(\ell_{\bar{i}2}-\ell_{\bar{i}3})+\ell_{i3}(2\ell_{\bar{i}3}-\ell_{\bar{i}2}-\ell_{\bar{i}4})\\\
0&0&\ell_{i5}-\ell_{i6}+\frac{\ell_{i6}(\ell_{i3}+\ell_{i6}-\ell_{i4}-\ell_{i5})}{1-\ell_{i4}+\ell_{i6}}&1+\ell_{i5}-\ell_{i3}\end{pmatrix}$
(55)

The solution of this smaller linear system reads

$\displaystyle\braket{\hat{n}_{\bar{i}\sigma}\hat{n}_{\bar{i}\bar{\sigma}}}=\tau_{i,1}\braket{\hat{n}_{\bar{i}\sigma}}$
(56a)
$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}}=\tau_{i,2}\braket{\hat{n}_{\bar{i}\sigma}}$
(56b)
$\displaystyle\braket{\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}=\tau_{i,3}\braket{\hat{n}_{\bar{i}\sigma}}$
(56c)
$\displaystyle\braket{\hat{n}_{i\sigma}\hat{n}_{i\bar{\sigma}}\hat{n}_{\bar{i}\sigma}}=\tau_{i,4}\braket{\hat{n}_{\bar{i}\sigma}}$
(56d)

where

$\displaystyle\tau_{i,1}=$
$\displaystyle\frac{1}{\tau_{i,0}}\left[\ell_{\bar{i}2}\left(\ell_{i3}-1\right)\left(\ell_{i5}-1\right)+2\left(1-\ell_{i3}\right)\ell_{\bar{i}3}\ell_{i5}+\ell_{i3}\ell_{\bar{i}4}\ell_{i5}\right]$
(57a) $\displaystyle\tau_{i,2}=$
$\displaystyle\frac{1}{\tau_{i,0}}\left[\ell_{\bar{i}2}\left(\ell_{i6}+\ell_{i5}\left(\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}-1\right)-\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\left(\ell_{i3}-\ell_{i4}+\ell_{i6}\right)\right)\right]$
$\displaystyle+\frac{1}{\tau_{i,0}}\left[1-\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\left(\ell_{\bar{i}2}+\ell_{\bar{i}4}-2\ell_{\bar{i}3}\right)\left(\ell_{i4}-\ell_{i3}\right)\right]\ell_{i5}$
(57b) $\displaystyle\tau_{i,3}=$
$\displaystyle\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\tau_{i,1}$ (57c)
$\displaystyle\tau_{i,4}=$
$\displaystyle\left(\ell_{i4}-\ell_{i3}\right)\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\tau_{i,1}+\ell_{i3}\tau_{i,2}$
(57d) $\displaystyle\tau_{i,0}=$
$\displaystyle\left(\ell_{\bar{i}2}(\ell_{i3}-2)+2(1-\ell_{i3})\ell_{\bar{i}3}+\ell_{i3}\ell_{\bar{i}4}\right)\left(\ell_{i5}-\ell_{i6}+\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\left(\ell_{i3}-\ell_{i4}-\ell_{i5}+\ell_{i6}\right)\right)$
$\displaystyle+\left[1-\frac{\ell_{i6}}{1-\ell_{i4}+\ell_{i6}}\left(\ell_{\bar{i}2}+\ell_{\bar{i}4}-2\ell_{\bar{i}3}\right)\left(\ell_{i4}-\ell_{i3}\right)\right](1+\ell_{i5}-\ell_{i3})$
(57e)

where the $\ell_{ij}$ are defined in Eq. 52a. Since in this way, all the
density correlators (and thus all the residues of the single-particle GF, see
Eq. 23) are expressed in terms of the densities themselves, the final step to
obtain the densities requires the solution of the $2\times 2$ linear system

$\displaystyle\begin{pmatrix}\eta_{11}&\eta_{12}\\\
\eta_{21}&\eta_{22}\end{pmatrix}\begin{pmatrix}\braket{\hat{n}_{i\sigma}}\\\
\braket{\hat{n}_{\bar{i}\sigma}}\end{pmatrix}=\begin{pmatrix}\phi(p_{1}^{1})\\\
\phi(p_{1}^{2})\end{pmatrix}$ (58)

where the parameters $\eta_{i,j}$ are given in Eq. 25. The solution of the
linear system finally gives the expression for the densities Eq. 24.

###### Acknowledgements.

We acknowledge useful discussions with Gianluca Stefanucci. We acknowledge
financial support through Grant PID2020-112811GB-I00 funded by
MCIN/AEI/10.13039/501100011033 as well as by grant IT1453-22 “Grupos
Consolidados UPV/EHU del Gobierno Vasco”. We acknowledge the technical support
provided by SGIker (Scientific Computing Services UPV/EHU).

## References

  * Hanson  _et al._ [2007] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Kloeffel and Loss [2013] C. Kloeffel and D. Loss, Annu. Rev. Condens. Matter Phys. 4, 51 (2013).
  * Loss and DiVincenzo [1998] D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 (1998).
  * Michler  _et al._ [2000] P. Michler, A. Kiraz, C. Becher, W. V. Schoenfeld, P. M. Petroff, L. Zhang, E. Y. Hu, and A. Imamoglu, Science 290, 2282 (2000).
  * Petta  _et al._ [2005] J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Science 309, 2180 (2005).
  * Reimann and Manninen [2002] S. M. Reimann and M. Manninen, Rev. Mod. Phys. 74, 1283 (2002).
  * Žutić  _et al._ [2004] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).
  * Veldhorst  _et al._ [2015] M. Veldhorst, C. H. Yang, B. Hwang, B. Huang, C. H. Yang, J. T. Muhonen, S. Simmons, A. Laucht, F. E. Hudson, R. Rahman, G. Klimeck, K. M. Itoh, A. Morello, and A. S. Dzurak, Nature 526, 410 (2015).
  * Zwanenburg  _et al._ [2013] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Rev. Mod. Phys. 85, 961 (2013).
  * Petta and Ralph [2001] J. R. Petta and D. C. Ralph, Phys. Rev. Lett. 87, 266801 (2001).
  * García de Arquer  _et al._ [2021] F. P. García de Arquer, D. V. Dmitri V. Talapin, V. I. Klimov, Y. Arakawa, M. Bayer, and E. H. Sargent, Science 373, eaaz8541 (2021).
  * Nygård  _et al._ [2000] J. Nygård, D. H. Cobden, and P. E. Lindelof, Nature 408, 342 (2000).
  * Buitelaar  _et al._ [2002] M. R. Buitelaar, A. Bachtold, T. Nussbaumer, M. Iqbal, and C. Schönenberger, Phys. Rev. Lett. 88, 156801 (2002).
  * Parks  _et al._ [2007] J. J. Parks, A. R. Champagne, G. R. Hutchison, S. Flores-Torres, H. D. Abruña, and D. C. Ralph, Phys. Rev. Lett. 99, 026601 (2007).
  * Anderson [1961] P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Hewson [1997] A. C. Hewson,  _The Kondo Problem to Heavy Fermions_ (Cambridge University Press, Cambridge, 1997).
  * Ono  _et al._ [2002] K. Ono, D. Austing, Y. Tokura, and S. Tarucha, Science 297, 1313 (2002).
  * van der Wiel  _et al._ [2002] W. G. van der Wiel, S. De Franceschi, J. Elzerman, T. Fujisawa, S. Tarucha, and L. Kouwenhoven, Rev. Mod. Phys. 75, 1 (2002).
  * You and Zheng [1999a] J. You and H. Zheng, Phys. Rev. B 60, 13314 (1999a).
  * Zhang  _et al._ [2023] Y. Zhang, S. Wang, W. Li, and M. Yue, Phys. Scr. 98, 105245 (2023).
  * Petta  _et al._ [2004] J. Petta, A. Johnson, C. Marcus, M. Hanson, and A. Gossard, Phys. Rev. Lett. 93, 186802 (2004).
  * Chen  _et al._ [2004] J. Chen, A. Chang, and M. Melloch, Phys. Rev. Lett. 92, 176801 (2004).
  * Fernández-Fernández  _et al._ [2022] D. Fernández-Fernández, Y. Ban, and G. Platero, Phys. Rev. Appl. 18, 054090 (2022).
  * Bouman  _et al._ [2020] D. Bouman, R. J. Van Gulik, G. Steffensen, D. Pataki, P. Boross, P. Krogstrup, J. Nygård, J. Paaske, A. Pályi, and A. Geresdi, Phys. Rev. B 102, 220505 (2020).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * Schollwöck [2011] U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Coleman [1984] P. Coleman, Phys. Rev. B 29, 3035 (1984).
  * Bickers [1987] N. E. Bickers, Rev. Mod. Phys. 59, 845 (1987).
  * Eckstein  _et al._ [2010] M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. 103, 056403 (2010).
  * Tanimura [2006] Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006).
  * Ishizaki and Fleming [2009] A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009).
  * Zubarev [1960] D. N. Zubarev, Sov. Phys. Usp. 3, 320 (1960).
  * Hubbard [1963] J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963).
  * Hubbard [1964a] J. Hubbard, Proc. R. Soc. Lond. A 277, 237 (1964a).
  * Hubbard [1964b] J. Hubbard, Proc. R. Soc. Lond. A 281, 401 (1964b).
  * Meir and Wingreen [1992] Y. Meir and N. S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992).
  * Van Roermund  _et al._ [2010] R. Van Roermund, S.-y. Shiau, and M. Lavagna, Phys. Rev. B 81, 165115 (2010).
  * Świrkowicz  _et al._ [2003] R. Świrkowicz, J. Barnaś, and M. Wilczyński, Phys. Rev. B 68, 195318 (2003).
  * You and Zheng [1999b] J. You and H.-Z. Zheng, Phys. Rev. B 60, 8727 (1999b).
  * Sztenkiel and Świrkowicz [2007] D. Sztenkiel and R. Świrkowicz, Acta Phys. Pol. A 111, 361 (2007).
  * Kuo and Chang [2007] D. M.-T. Kuo and Y.-C. Chang, Phys. Rev. Lett. 99, 086803 (2007).
  * Lamba and Joshi [2000] S. Lamba and S. Joshi, Phys. Rev. B 62, 1580 (2000).
  * Palacios  _et al._ [1997] J. Palacios, L. Liu, and D. Yoshioka, Phys. Rev. B 55, 15735 (1997).
  * Sierra  _et al._ [2016] M. A. Sierra, M. Saiz-Bretín, F. Domínguez-Adame, and D. Sánchez, Phys. Rev. B 93, 235452 (2016).
  * Bułka and Kostyrko [2004] B. R. Bułka and T. Kostyrko, Phys. Rev. B 70, 205333 (2004).
  * Sun and Guo [2002] Q.-F. Sun and H. Guo, Phys. Rev. B 66, 155308 (2002).
  * Chi and Li [2006] F. Chi and S.-S. Li, J. Appl. Phys. 99, 043705 (2006).
  * Chang and Kuo [2008] Y.-C. Chang and D. M.-T. Kuo, Phys. Rev. B 77, 245412 (2008).
  * Feng and Oppeneer [2011] Q. Feng and P. M. Oppeneer, J. Phys.: Condens. Matter 23, 425601 (2011).
  * Hewson [1966] A. Hewson, Phys. Rev. 144, 420 (1966).
  * Ryndyk  _et al._ [2009] D. A. Ryndyk, R. Gutiérrez, B. Song, and G. Cuniberti, Green function techniques in the treatment of quantum transport at the molecular scale, in  _Energy Transfer Dynamics in Biomaterial Systems_ (Springer, Berlin, 2009) pp. 213–335.
  * Haug and Jauho [2008] H. Haug and A.-P. Jauho,  _Quantum kinetics in transport and optics of semiconductors_ , Vol. 2 (Springer, Berlin, 2008).
  * Niu  _et al._ [1999] C. Niu, D. Lin, and T. Lin, J. Phys.: Condens. Matter 11, 1511 (1999).
  * Stefanucci and Van Leeuwen [2013] G. Stefanucci and R. Van Leeuwen,  _Nonequilibrium many-body theory of quantum systems: a modern introduction_ (Cambridge University Press, Cambridge, 2013).
  * Beenakker [1991] C. W. J. Beenakker, Phys. Rev. B 44, 1646 (1991).
  * Liul  _et al._ [2023] M. Liul, A. Ryzhov, and S. Shevchenko, Eur. Phys. J. Spec. Top. 232, 3227 (2023).
  * Mazal  _et al._ [2019] Y. Mazal, Y. Meir, and Y. Dubi, Phys. Rev. B 99, 075433 (2019).
  * Sobrino  _et al._ [2021] N. Sobrino, F. Eich, G. Stefanucci, R. D’Agosta, and S. Kurth, Phys. Rev. B 104, 125115 (2021).
  * Sobrino  _et al._ [2020] N. Sobrino, S. Kurth, and D. Jacob, Phys. Rev. B 102, 035159 (2020).
  * Tanimura [1990] Y. Tanimura, Phys. Rev. A 41, 6676 (1990).
  * Yan  _et al._ [2004] Y.-A. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395, 216 (2004).
  * Xu  _et al._ [2005] R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo, and Y. Yan, J. Chem. Phys. 122, 041103 (2005).
  * Sobrino [2024] N. Sobrino, EOM_DQD, https://github.com/Nahualcsc/EOM_DQD (2024).
  * Stefanucci and Kurth [2015] G. Stefanucci and S. Kurth, Nano Lett. 15, 8020 (2015).
  * Kurth and Stefanucci [2017] S. Kurth and G. Stefanucci, J. Phys.: Condens. Matter 29, 413002 (2017).
  * Jacob and Kurth [2018] D. Jacob and S. Kurth, Nano Lett. 18, 2086 (2018).
  * Sobrino  _et al._ [2019] N. Sobrino, R. D’Agosta, and S. Kurth, Phys. Rev. B 100, 195142 (2019).
  * Sobrino  _et al._ [2023] N. Sobrino, R. D’Agosta, and S. Kurth, Phys. Rev. B 107, 195116 (2023).