*3.1. Transient Charging of Excited States*

We consider a two-dimensional system of length *Lx* and width *Ly* described by a lattice with *Nx* sites on the *x* axis and *Ny* sites on the *y* axis. The total number of sites is denoted by *Nxy* = *NxNy*. By setting the two lattice constants *ax* and *ay* one has *Lx* = *axNx* and *Ly* = *ayNy*. Once we know the single-particle eigenstates of the electronic subsytem *S*1 we can write down its Hamiltonian *HS*1 := *H*(0) *S*1 + *W* in a second quantized form w.r.t. this basis, that is:

$$\begin{array}{rcl} H\_{\mathbb{S}\_1}^{(0)} & = & \sum\_{n,n'=1}^{N\_{\text{xy}}} \langle \psi\_{n} | \hat{h}\_{\text{S}\_1}^{(0)} | \psi\_{n'} \rangle c\_{n}^{\dagger} c\_{n'} = \sum\_{n} \varepsilon\_{n} c\_{n}^{\dagger} c\_{n} \\ \mathcal{W} & = & \frac{1}{2} \sum\_{n,m,n',m'} V\_{nmn'm'} \mathfrak{c}\_{n}^{\dagger} \mathfrak{c}\_{m}^{\dagger} \mathfrak{c}\_{m'} \mathfrak{c}\_{n'} \end{array} \tag{41}$$

where the Coulomb matrix elements are given by (**<sup>r</sup>**,**r** are sites of the 2D lattice):

$$V\_{nmn'm'} = \sum\_{\mathbf{r},\mathbf{r}'} \psi\_n^\*(\mathbf{r}) \psi\_{m}^\*(\mathbf{r}') V\_{\mathbb{C}}(\mathbf{r} - \mathbf{r}') \psi\_{n'}(\mathbf{r}) \psi\_{m'}(\mathbf{r}'). \tag{42}$$

The Coulomb potential itself is given by

$$V\_{\mathbb{C}}(\mathbf{r}, \mathbf{r}') = \frac{e^2}{4\pi\epsilon(|\mathbf{r} - \mathbf{r}'| + \eta)}\,,\tag{43}$$

where *η* is a small positive regularization parameter.

Like in the continuous model, the tunneling coefficients *Tlqn* are associated to a pair of states {*ψlq*, *ψn*} from the lead *l* and the sample *S*1. However the lattice version is much simpler:

$$T\_{q\mathbf{u}}^{l} = V\_{l} \psi\_{q}^{l\*} (0\_{l}) \psi\_{\mathbf{u}}(i\_{l}) , \tag{44}$$

where 0*l* is the site of the lead *l* which couples to the contact site *il* in the sample. The wavefunctions of the semi-infinite lead are known analytically:

$$\psi\_q^l(j) = \frac{\sin(q(j+1))}{\sqrt{2\pi \sin q}}, \quad \varepsilon\_q = 2\pi \cos q. \tag{45}$$

In the above equation *τ* is the hopping energy of the leads. The integral over *q* in the tunneling Hamiltonian (see Equation (8) from Section 2) counts the momenta of the incident electrons such that *ε<sup>l</sup>*(*q*) scans the continuous spectrum of the semi-infinite leads *σl* ∈ [−2*<sup>τ</sup>* + Δ, 2*τ* + Δ] where Δ is a shift which is chosen such that *σl* covers the lowest-energy many-body spectrum of the central system. The construction of the coupling coefficients *Tlqn* shows that a single-particle state which vanishes at the contact sites does not contribute to the currents. This is the case for states which are mostly localized at the center of the sample, while in the presence of a strong magnetic field the currents will be carried by edge states.

In [12] we implemented GME for a non-interacting lattice Hamiltonian, whereas the Coulomb interaction effects were introduced in [14]. In what concerns the geometrical effects we essentially showed that the transient currents depend on the location of the contacts (through the value of the single-particle wavefunctions of the sample at those points) but also on the initial state and on the switching functions *<sup>χ</sup>l*(*t*) of the leads. It turns out that the stationary current does not depend on the last two parameters, in agreemen<sup>t</sup> with rigorous results [60,61]. We also identified a delay of the output currents which was attributed to the electronic propagation time along the edge states of the Hofstadter spectrum.

The presence of Coulomb interaction brings in specific steady-state features known from previous calculations like the Coulomb blockade and the step-like structure of the current–voltage characteristics. On the other hand the GME method naturally allows a detailed analysis of the time-dependent currents associated to each many-body configuration as well as of the relevant populations.

Since the Hamiltonian *HS*1 of the interacting system commutes with the total number operator *Q* = ∑*n <sup>c</sup>*†*ncn*, its eigenstates |*ν* can still be labelled by the occupation *Nν* of the non-interacting MBSs from which the state is built. Then the single index *ν* can be replaced by two indices, the particle number *Nν* and an index *iν* = 0, 1, 2, ... for the ground (*iν* = 0) and excited states (*iν* > 0, where *iν* also counts the spin degeneracy). The notation for the interacting many-body energies is changed accordingly E*ν* → E(*<sup>i</sup>ν*) *Nν*.

We now define some useful quantities for our time-dependent analysis. The charge accumulated on *N*-electrons states is calculated by collecting the associated populations:

$$q\_N(t) = \varepsilon N \sum\_{\nu, \mathbf{n}\_y = N} \langle \nu | \rho(t) | \nu \rangle \,\tag{46}$$

where the sum counts all states whose total occupation *nν* = *N*. Similarly one can identify the transient currents *Jl*,*<sup>N</sup>* carried by *<sup>N</sup>*-particle states. These currents can be traced back form the right hand side of the GME:

$$
\langle \dot{Q} \rangle = \sum\_{N} (f\_{\mathbb{L},N}(t) - f\_{\mathbb{R},N}(t)) = \sum\_{N} \dot{q}\_{N}(t). \tag{47}
$$

Throughout this work we shall adopt the following sign convention for the currents associated to each lead: *JL* > 0 if the electrons flow from the left lead towards the sample and *JR* > 0 if they flow from the sample towards the right lead.

The sequential tunneling processes change the many-body configurations of the electronic system. The energy required to bring the system to the *i*-th MBS with *N* particles is measured w.r.t. the ground state with *N* − 1 electrons (*i* = 0, 1, ...). We introduce two classes of chemical potentials of the sample:

$$
\mu\_{\mathcal{g},N}^{(i)} = \mathcal{E}\_N^{(i)} - \mathcal{E}\_{N-1}^{(0)} \tag{48}
$$

$$
\mu\_{x,N}^{(i)} = \mathcal{E}\_N^{(i)} - \mathcal{E}\_{N-1}^{(1)} \tag{49}
$$

where *μ*(*i*) *g*,*N* characterizes transitions from the ground state (*N* − <sup>1</sup>)-particle configuration to various *<sup>N</sup>*-particle configurations. In particular *μ*(0) *g*,*N* describes addition processes involving ground-states with *N* − 1 and *N* electrons while *μ*(*<sup>i</sup>*><sup>0</sup>) *g*,*N* refers to transitions from (*N* − <sup>1</sup>)-particle ground state to excited *<sup>N</sup>*-particle configurations. The chemical potentials *μ*(*i*) *<sup>x</sup>*,*N* describe transitions from the 1st (*N* − <sup>1</sup>)-particle excited states to configurations with *N* particles. In a transition of this type an electron tunnels on the lowest single particle state to the central system which already contains one electron on the excited single-particle state |*<sup>σ</sup>*2. As a result some of the triplet states are being populated. We shall see that these transitions play a role especially in the transient regime.

For numerical calculations we considered a 2D quantum wire of lenght *Lx* = 75 nm and width *Ly* = 10 nm. The lowest two spin-degenerate single-particle levels are *ε*1 = 0.375 meV and *ε*2 = 3.37 meV. The non-interacting MBSs are described by the spins of the occupied single-particle levels, e.g., |*<sup>σ</sup>*1*σ*2 is a two-particle configuration with a spin *σ* associated to the lowest single-particle state and a second electron with opposite spin orientation on the energy level *ε*2. Besides the usual singlet (*S*) and triplet (*T*) states we find that the Coulomb interaction induces the configuration mixing of the antiparallel configurations | ↑1↓1 and | ↑2↓2. More precisely, we ge<sup>t</sup> an interacting ground two-particle state mostly made of | ↑1↓1 (whose weight is 0.86) and with a small component (weight 0.14) of state | ↑2↓2. Conversely, | ↑1↓1 is also found in the highest energy two-particle state. We stress here that the configuration mixing decreases and eventually vanishes if the gap *<sup>E</sup>*↑2↓2 − *<sup>E</sup>*↑1↓1 between two non-interacting energies is much larger that the corresponding matrix element.

In Figure 2a we show the chemical potentials corresponding to interacting MB configurations with up to three electrons. As long as the chemical potential *μ*(*i*) *g*,*N* lies within the bias window the corresponding state will contribute both to the transient and steady-state currents. We shall see that if *μ*(*i*) *g*,*N* < *μR* the state |*i*, *N* contributes only to the transient currents. Finally, when *μ*(*i*) *g*,*N μL* the state |*i*, *N* is poorly populated and will not contribute to transport. Let us stress here a rather unusual transition from |*<sup>σ</sup>*2 to the ground two-particle state which is mostly made of |*<sup>σ</sup>*1*σ*1. The corresponding addition energy *μ*(0) *<sup>x</sup>*,2 = 2 meV is smaller than the energy required for the usual transition |*<sup>σ</sup>*1→|*<sup>σ</sup>*1*σ*1. This happens because of the Coulomb mixing between | ↑1↓1 and | ↑2↓2 which makes possible the transition from the excited single-particle state to the mixed interacting two-particle groundstate. The chemical potential *μ*(2) *<sup>x</sup>*,2 describes the transition from the excited single-particle state |*<sup>σ</sup>*2 to the triplet states.

**Figure 2.** (**a**) The chemical potentials *μ*(*i*) *g*,*N* (red crosses) and *μ*(*i*) *<sup>x</sup>*,*N* (blue crosses) for *<sup>N</sup>*-particle configurations, *N* = 1, 2, 3. For a given particle number *N* the chemical potentials are ordered vertically according to the index *i* = 0, 1, .... The horizontal lines correspond to specific values of the chemical potentials in the leads (see the discussion in the text); (**b**) The time-dependent currents in the left lead at different values of the bias window *μL* − *μR*.

Already by analysing Figure 2 one can anticipate to some extend how the transport takes place in terms of allowed sequential tunneling processes. Suppose that the chemical potentials of the leads are selected such that *μ*(0) *g*,1 < *μR* < *μ*(1) *g*,1 < *μL* < *μ*(0) *g*,2 (as an example we set *μR* = 1 meV and *μL* = 4 meV). Then both single-particle levels are available for tunneling but one expects that the double occupancy is excluded because *μL* < *μ*(0) *g*,2 ∼ 5 meV. According to this scenario, more charge will accumulate on *ε*1, the excited states |*<sup>σ</sup>*2 will eventually deplete and the steady-state current vanishes in the steady-state. This is the well known Coulomb blockade effect. However, we see in Figure 2b that the steady-state current vanish only when *μL* < *μ*(1) *g*,1 as well, which sugges<sup>t</sup> that the presence of the excited single-particle states within the bias window leads to a partial lifting of the Coulomb blockade. We stress that such an effect cannot be predicted within a single-site model with onsite Coulomb interaction. A third curve shows the current for *μL* = 5.5 meV and *μL* = 4 meV.

Figure 3 presents the evolution of the relevant populations at two values of the bias window. In Figure 3a the population *<sup>P</sup>*1*g* = *<sup>P</sup>*↑1 + *<sup>P</sup>*↓1 of the ground single-particle states dominates in the steady state. This is expected, as the corresponding chemical potential lies below the bias window so this state will be substantially populated. The other configurations contributing to the steady-state are just the ones which can be populated by tunnelings from the left lead, that is the excited single-particle states and all two-particle states except for the single configuration which cannot be accessed. By looking at Figure 2a one infers that the two-particle states are being populated when one more electron is added

from the left lead on the initial excited single-particle state |*<sup>σ</sup>*2. In particular, the ground two-particle state is populated only due to the Coulomb-induced configuration mixing.

**Figure 3.** The populations of ground (g) and excited (x) *<sup>N</sup>*-particle states (*N* = 1, 2, 3) for two bias windows; (**a**) *μL* = 4 meV, *μR* = 2 meV; (**b**) *μL* = 5.5 meV, *μR* = 4 meV. In panel (**a**) *P*3 is negligible and was omitted.

A completely different behavior is noticed in Figure 3b. As the bias window is pushed upwards such that *μ*(1) *g*,1 < *μR* < *μ*(1) *g*,2 < *μL* the transitions from the lowest states |*<sup>σ</sup>*1 to two-particle states are also activated. Consequently, the population *P*2*x* of the excited two-particle configurations exceeds *<sup>P</sup>*1*g* and dominate in the steady-state regime. Note that *P*2*x* > *<sup>P</sup>*2*g* because it collects the population of the degenerate triplet states. In the transient regime the excited single particle states are populated much faster than the ground states. This happens because of the different localizations of the single-particle wavefunctions on the contact regions. We find that the wavefunction associated to the 2nd single-particle state has a larger value at the endpoints of the leads. A drop of *P*1*x* follows as the ground one-electron states and the other two-particle configurations become active (a similar feature is noticed in Figure 3a). A small populations of the three particle states can be also observed. The steady-state current increases considerably (see Figure 2b) and is due to the two-particle states.

We end this section with a discussion on the partial currents *JL*,*<sup>N</sup>* and *JR*,*<sup>N</sup>* associated to *<sup>N</sup>*-particle states. Although they cannot be individually measured, these currents provide further insight into the transport processes, in particular on the way in which the steady-state regime is achieved.

Figure 4a shows that in the steady-state regime the currents carried by the one-particle states *JL*,<sup>1</sup> and *JR*,<sup>1</sup> achieve a negative value when *μ*(0) *g*,1 < *μ*(1) *g*,1 < *μR* < *μ*(0) *g*,2 < *μ<sup>L</sup>*, whereas the two-particle currents evolve to a larger positive value such that the total current *JL* will be positive as already shown in Figure 2b. When the bias window is shifted down to *μL* = 4 meV and *μR* = 2 meV all transients are mostly positive (see Figure 4b). One observes that the single-particle configurations are responsible for the spikes of the total current *JL* and that the two-particle currents display a smooth behavior. These features can be explained by looking at the charge occupations *qN* shown in Figure 4c,d. As *μ*(0) *g*,1 and *μ*(1) *g*,1 are both below *μR* = 4 meV while *μ*(0) *g*,2 is well within the bias window, the population of the single-particle states increases rapidly in the transient regime but then also drops in favour of *P*2, the total occupation of two-particle states. Such a redistribution of charge among configurations with different particle numbers is less pronounced in Figure 4d, because in this case the smaller contribution of the two-particle states is only due to transitions allowed by *μ*(0) *<sup>x</sup>*,2 and *μ*(1) *<sup>x</sup>*,2 which are now located within the bias window. The slope of *q*2 also changes sign in the transient regime and one can check from Figure 4b that on the corresponding time range *JR*,<sup>2</sup> slightly exceeds *JL*,2.

**Figure 4.** (**a**) The transient currents *JL*,*<sup>N</sup>* and *JR*,*<sup>N</sup>* associated to one and two-particle configurations for *μL* = 5.5 meV, *μR* = 4 meV; (**b**) The same currents for a bias window *μL* = 4 meV, *μR* = 2 meV; (**c**) The charge *qN* accumulated on *<sup>N</sup>*-particle states at *μL* = 5.5 meV, *μR* = 4 meV; (**d**) *qN* for *μL* = 4 meV, *μR* = 2 meV, and *qN* are given in units of electron charge *e*.

The occupation of the three-particle configurations is negligible so *q*3 is also small and was included here only for completeness while the associated currents were omitted.

### *3.2. Coulomb Switching of Transport in Parallel Quantum Dots*

After using the GME formalism to describe transient transport via excited states in a single interacting nanowire we now extend its applications to capacitively coupled quantum systems. Besides Coulomb blockade, the electron-electron interaction cause momentum-exchange which leads to the well known Coulomb drag effect in double-layer structures [62] and double quantum dots [63–65] or wires [66]. Also, theoretical calculations on thermal drag between Coulomb-coupled systems were recently presented [67,68].

Here we consider a very simple model for two parallel quantum dots [17] (a sketch of the system is given in Figure 5). Each system is described by a 1D four-sites chain and for simplicity we neglect the spin degree of freedom which will only complicate the discussion of the effects. The diagonalization procedure provides all 256 many-body configurations emerging from the 8 single-particle states. Let us point out that the interdot and intradot interactions are treated on equal footing beyond the single-capcitance model. The hopping energy within the dots is *tD* = 1 meV and the time unit is expressed in units of *h*¯ /*tD*. Then the currents are calculated in units of *etD*/¯*h*. The tunneling rates to the four leads are all equal *VLa* = *VRa* = *VLb* = *VRb*.

We shall use the GME method to study the onset of the interdot Coulomb interaction. In order to distinguish the transient features due to mutual capacitive coupling we consider a transport setting in which each dot is connected to the leads at different times. More precisely, one system, say *QDa* is open at the initial instant *ta* = 0 and then reaches a stationary state (*JLa* = *JRa*) at some later time *Ta*. The coupling of the nearby system to its leads is switched on at *tb* > *Ta* such that the changes in the current *Ja* can only be due to mutual Coulomb interaction. Note that the usual Markov–Lindblad

version of the master equation simulate the transport when the four leads are coupled suddenly and simultaneously to the double-dot structure.

**Figure 5.** A sketch of the parallel double-dot system. Each QD is coupled to source-drain particle reservoirs described by chemical potentials *μLs* and *μRs*, *s* = *a*, *b*. There is no interdot electron tunneling but the systems are correlated via Coulomb interaction.

As before, the interacting many-body configurations can be labeled by to the occupations of each dot according to the correspondence E*ν* → E(*<sup>i</sup>ν*) *Na*,*ν*,*Nb*,*<sup>ν</sup>* . Here *Ns*,*<sup>ν</sup>* is the number of electrons in the system *s* associated to a many-body configuration *ν*. If the two systems are identical the lowest chemical potentials are introduced as:

$$
\mu^{(0)}\_{\mathcal{S}}(N\_{\mathfrak{a}}, N\_{\mathfrak{b}}) = \mathcal{E}^{(0)}\_{N\_{\mathfrak{a}}, N\_{\mathfrak{b}}} - \mathcal{E}^{(0)}\_{N\_{\mathfrak{a}}-1, N\_{\mathfrak{b}}} = \mathcal{E}^{(0)}\_{N\_{\mathfrak{a}}, N\_{\mathfrak{b}}} - \mathcal{E}^{(0)}\_{N\_{\mathfrak{a}}, N\_{\mathfrak{b}}-1} \tag{50}
$$

.

because of the degeneracy w.r.t. to the total occupation number E(0) *Na*,*Nb* = E(0) *Nb*,*Na* For the parameters chosen here one finds: *μ*(0) *g* (1, 1) = 3 meV, *μ*(0) *g* (2, 0) = 4 meV and *μ*(0) *g* (2, 1) = 4.5 meV. The location of the several chemical potentials w.r.t. the two bias windows already suggests the possible interdot correlation effects. The main point is that the transport channels through one dot also depend on the occupation of the nearby dot. One therefore expects that the currents *JLa* and *JRa* also depend on the bias applied on the nearby system.

To discuss this effect we performed transport simulations for two arrangements (A and B) of the bias window *μLs* − *μRs*. In the A-setup we select the four chemical potentials such that the chemical potentials associated to the many-body configurations relevant for transport obey the inequalities *μ*(0) *g* (1, 1) < *μRs* < *μ*(0) *g* (2, 0) < *μLs* < *μ*(0) *g* (2, 1). The scenario is easy to grasp: As *QDa* is coupled to the leads and the nearby dot is disconnected and empty, it will accumulate charge and evolve to a steady state where the current is essentially given by tunneling assisted transitions between E(0) 2,0 ↔ E(0) 1,0 This behavior is observed in Figure 6a up to *tb* = 150 ps when *QDb* is also coupled to its leads. Note also that the charge occupation of *QDa* almost saturates at *Qa* = 1.6. As expected, for *t* > *tb* a transient current develops in *QDb*, but a simultaneous drop the *JLa* and *JRa* shows the dynamical onset of the charge sensing effect between the two systems. In the final steady-state the two currents nearly vanish, thus proves their negative correlation due to the mutual Coulomb interaction. The charges *Qa*,*<sup>b</sup>* reach the same value and sugges<sup>t</sup> that in the long time limit the double system contains one electron on each dot. Remark that in the final steady-state the dominant population corresponds to the many-body energy E(0) 1,1 which is not favorable for transport through any of the dots as long as *μ*(0) *g* (2, 1) = E(0) 2,1− E(0) 1,1is outside the bias window.

Figure 6b,d present the currents and the charge occupations for the second setup B which is defined by the inequalities *μ*(0) *g* (2, 0) < *μRs* < *μ*(0) *g* (2, 1) < *μLs*. Following the same reasoning as before one infers that now *QDa* will enter the Coulomb blockade regime before *t* = *tb* because there are no transport channels within the bias window. However, the blockade is removed due to the second dot whose charging activates tunneling through *μ*(0) *g* (2, 1) = *μ*(0) *g* (1, <sup>2</sup>). This is an example of positive correlations between the two systems. Further discussions can be found in a previous publication [17].

**Figure 6.** The transient currents in the two systems for different chemical potentials of the leads: (**a**) *μLa* = *μLb* = 4.25 meV, *μRa* = *μRb* = 3.75 meV; (**b**) *μLa* = *μLb* = 4.75 meV, *μRa* = *μRb* = 4.35 meV; (**c,d**) The charge occupations of the two systems associated to the currents in Figure 6a,b. The charges *Qa*,*<sup>b</sup>* are given in units of electron charge *e*.
