*2.3. Numerical Implementation and Observables*

The last step before numerical implementation requires the calculation of the system-environment couplings *HT* and *HE* w.r.t. the full basis |*ϕp*). Clearly, to this end we shall use the unitary transformations |*λ*↔|*ν* and |*ν*, *j*↔|*ϕp*) which are already at hand due to the stepwise

diagonalization procedure introduced in the previous section. Then let us introduce some generalized 'jump' operators collecting all transitions, between pairs of fully interacting states, generated by tunneling of an electron with momentum *q* from the *l*-th lead to the single-particle levels of the electronic system *S*1:

$$\mathcal{T}\_l(q) = \sum\_{p, p'} \mathcal{T}\_{pp'}^{l}(q) |\varphi\_p\rangle\langle\varphi\_{p'}| \; , \quad (\mathcal{T}\_l(q))\_{pp'} = \sum\_n T\_{qn}^{l}(\varphi\_p|c\_n^\dagger|\varphi\_{p'}) \; . \tag{30}$$

Then the dissipation operator associated to the particle reservoirs reads:

$$\mathcal{D}\_{\text{leads}}[\rho, t] = -\frac{1}{\hbar^2} \sum\_{l=\text{L,R}} \int dq \,\chi\_l(t) \left( [\mathcal{T}\_l, \Omega\_{\text{ql}}(t)] + h.c.\right), \tag{31}$$

with the following notation:

$$
\Omega\_{ql}(t) = \mathcal{U}\_{\mathcal{S}}(t, t\_0) \int\_{t\_0}^{t} ds \, \chi\_{l}(s) \Pi\_{ql}(s) e^{i(\left(s - t\right)/\hbar\right) \mathbf{c}\_{l}(q)} \mathcal{U}\_{\mathcal{S}}^{\dagger}(t, t\_0), \tag{32}
$$

$$\Pi\_{ql}(\mathbf{s}) = \mathcal{U}\_{\mathbf{S}}^{\dagger}(\mathbf{s}, t\_0) \left( \mathcal{T}\_l^{\dagger} \rho(\mathbf{s}) (1 - f\_l) - \rho(\mathbf{s}) \mathcal{T}\_l^{\dagger} f\_l \right) \mathcal{U}\_{\mathbf{S}}(\mathbf{s}, t\_0), \tag{33}$$

and where for simplicity we omit to write the energy dependence of the Fermi function *fl*. Similarly, the bosonic operators have to written down w.r.t. the full basis which then leads to the calculation of Dbath. Under the Markov approximation w.r.t. the correlation function of the bosonic reservoir the latter becomes local in time.

The GME is solved numerically by time discretization using the Crank–Nicholson method which allows us to compute the reduced density operator for discrete time steps *ρ*(*tn*), starting with an initial condition corresponding to a given state of the isolated central system, i.e., before the onset of the coupling with the leads. We take advantage of the fact that, by discretizing the time domain, the operator <sup>Ω</sup>*q<sup>l</sup>*(*tn*+<sup>1</sup>) obeys a recursive formula generated by the incremental integration between *tn* and *tn*+1, that is:

$$
\Omega\_{ql}(t\_{n+1}) = l L\_{\mathbb{S}}(t\_{n+1}, t\_n) \Omega\_{ql}(t\_n) l L\_{\mathbb{S}}^{\uparrow}(t\_{n+1}, t\_n) + \mathcal{A}\_{ql}(t\_{n+1}, t\_n; \rho(t\_{n+1}), \rho(t\_n)), \tag{34}
$$

where the second term of the right-hand side depends on the ye<sup>t</sup> unknown *ρ*(*tn*+<sup>1</sup>). For any pair of time steps {*tn*, *tn*+<sup>1</sup>} we initially approximate *ρ*(*tn*+<sup>1</sup>) in <sup>A</sup>*q<sup>l</sup>* by the already calculated *ρ*(*tn*), and perform iterations to recalculate *ρ*(*tn*+<sup>1</sup>) via the GME, each time updating *ρ*(*tn*+<sup>1</sup>) in <sup>A</sup>*ql*, until a convergence test for *ρ*(*tn*+<sup>1</sup>) is fulfilled. At any step of the iteration we also calculate and include <sup>D</sup>bath[*ρ*, *tm*] into the iterative procedure; its calculation is much simpler as the Markov approximation w.r.t. the bath degrees of freedom takes care of the time integral so this dissipative term becomes local in time. Finally, we check numerically the conservation of probability and the positivity of the diagonal elements of *ρ*(*tm*), i.e., the populations of fully interacting states |*ϕp*) at the corresponding time step and for each iteration.

There are several reasons to extend the GME method beyond single-level models. (1) The electronic transport at finite bias collects contributions from all the levels within the bias window. This feature leads to the well known stepwise structure of the current-voltage characteristics; (2) In the presence of Coulomb interaction the GME must be derived in the language of many-body states which allows us to perform exact diagonalization on appropriate Fock subspaces; (3) The minimal model which describes the effect of the field-matter coupling in optical cavities with embedded quantum dots requires at least two single-particle levels.

Both the GME and non-equilibrium Green's function formalism (NEGF) rely on the partitioning approach and allow for many-body interaction in the central system, while the leads are assumed to be non-interacting (this assumption leads in particular to the Fermi distribution of the particle reservoirs). There is however a crucial difference between the two methods. The perturbative expansion of the

dissipative kernel forces restricts the master equation approach to weak lead-sample tunnelings while the interaction effects are accounted for exactly. In contrast, the Keldysh formalism is not limited to small system-reservoir couplings but the Coulomb effects have to be calculated from appropriate interaction self-energies. Which method fits better is simply decided by the particular problem at hand.

As stated in the Introduction, the advantage of the RDO stems from the fact that it can be used to calculate statistical averages of various observables O of the hybrid system:

$$
\langle \mathcal{O} \rangle = \text{Tr}\_{\mathcal{F}} \{ \rho(t) \mathcal{O} \} \,. \tag{35}
$$

Useful examples are averages of the photon number operator N *ph* = *a*†*a* and of the charge operator Q = ∑*n <sup>c</sup>*†*ncn*. Also, the average currents in a two-lead geometry (i.e., *l* = *L*, *R* can be identified from the continuity equation:

$$
\langle Q \rangle = \text{Tr}\_{\mathcal{F}} \{ \mathcal{Q}\phi(t) \} = \mathcal{J}\_L(t) - \mathcal{J}\_R(t). \tag{36}
$$

### *2.4. Coupling between Leads and Central System*

The modeling of the central systems and the reservoirs can be performed either by using continuous confining potentials or a spatial grid. Examples are a short parabolic wire [44,45], ring [46,47], parallel wires with a window coupler [48], and wire with embedded dot [44,49] or dots [50]. The coupling between the leads and the central system with length *Lx* is described by Equation (8), and in order to reproduce scattering effects seen in a Lippmann–Schwinger formalism [15,51,52] the coupling tensor is defined as

$$T\_{qn}^{l} = \int\_{\Omega\_{\mathcal{S}}^{l} \times \Omega\_{l}} d\mathbf{r} d\mathbf{r}' \left(\Psi\_{q}^{l}(\mathbf{r}')\right)^{\*} \Psi\_{n}^{S}(\mathbf{r}) \mathbf{g}\_{qn}^{l}(\mathbf{r}, \mathbf{r}') + h.c.,\tag{37}$$

for states with wavefunction Ψ*lq* in lead *l*, and Ψ*Sn* in the central system. The domains for the integration of the wavefunctions in the leads are chosen to be

$$\begin{array}{rcl}\Omega\_{\text{L}} &=& \left\{ (\mathbf{x}, \mathbf{y}) \, \middle| \, \left[ -\frac{L\_{\mathbf{r}}}{2} - 2a\_{\mathbf{w}\prime} - \frac{L\_{\mathbf{r}}}{2} \right] \times \left[ -3a\_{\mathbf{w}\prime} + 3a\_{\mathbf{w}\prime} \right] \right\}, \\ \Omega\_{\text{R}} &=& \left\{ (\mathbf{x}, \mathbf{y}) \, \middle| \, \left[ +\frac{L\_{\mathbf{r}}}{2}, +\frac{L\_{\mathbf{r}}}{2} + 2a\_{\mathbf{w}\prime} \right] \times \left[ -3a\_{\mathbf{w}\prime} + 3a\_{\mathbf{w}\prime} \right] \right\}, \end{array} \tag{38}$$

and for the system as

$$\begin{array}{rcl}\Omega\_{\mathbb{S}}^{\mathrm{L}} &=& \left\{ (\mathbf{x}, \mathbf{y}) \,|\, \left[ -\frac{L\_{\mathbf{x}}}{2}, -\frac{L\_{\mathbf{x}}}{2} + 2a\_{\mathrm{w}} \right] \times \left[ -3a\_{\mathrm{w}\prime} + 3a\_{\mathrm{w}} \right] \right\},\\\Omega\_{\mathbb{S}}^{\mathrm{R}} &=& \left\{ (\mathbf{x}, \mathbf{y}) \,|\, \left[ +\frac{L\_{\mathbf{x}}}{2} - 2a\_{\mathrm{w}\prime} + \frac{L\_{\mathbf{x}}}{2} \right] \times \left[ -3a\_{\mathrm{w}\prime} + 3a\_{\mathrm{w}} \right] \right\}. \end{array} \tag{39}$$

The function

$$\mathcal{g}\_{q\eta}^{l}(\mathbf{r},\mathbf{r}') = \mathcal{g}\_0^l \exp\left[-\delta\_1^l(\mathbf{x}-\mathbf{x}')^2 - \delta\_2^l(y-y')^2\right] \exp\left(\frac{-|E\_n - \epsilon^l(q)|}{\Delta\_E^l}\right). \tag{40}$$

with **r** ∈ Ω*l*S and **r** ∈ Ω*l* determines the coupling of any two single-electron states by the "nonlocal overlap" of their wave functions in the contact region of the leads and the system, and their energy affinity. A schematic view of the coupling is seen in Figure 1. The parameters *δl*1 and *δl*1 define the spatial range of the coupling within the domains Ω*l*S × Ω*l* [44].

**Figure 1.** A schematic of the coupling of the system to the leads. The transparent green areas correspond to the contact regions defined by the nonlocal overlap function *g*L,R *qn* in *HT*(*t*).

The short quantum wire is considered to have hard wall confinement in the transport direction, the *x*-direction, at *x* = ±*Lx*/2, and parabolic confinement in the *y*-direction with characteristic energy *h*¯ Ω. Possibly, the leads and the central system are considered to be placed in a perpendicular homogeneous external magnetic field **B** = *B***z**ˆ. Together they lead to a natural length scale, the effective magnetic length *aw* with *<sup>a</sup>*2*w*Ω*<sup>w</sup>* = *h*¯ /*<sup>m</sup>*, with <sup>Ω</sup>2*w* = [(<sup>Ω</sup>0)<sup>2</sup> + (*<sup>ω</sup>s*)<sup>2</sup>]1/2, and the cyclotron frequency *ωc* = (*eB*/*m*). For GaAs with effective mass *m* = 0.067*me* relative dielectric constant *κc* = 12.3 and confinement energy ¯*h*Ω0 = 2.0 meV, *aw* = 23.8 nm. The magnetic field *B* = 0.1 T.

The energy spectrum of the quasi 1D semi-infinite lead *l* is represented by  *<sup>l</sup>*(*q*), with *q* standing for the momentum quantum number of a standing wave state, and the subband index *<sup>n</sup>l*. The spectrum in the absence of spin orbit interactions can be evaluated exactly analytically [53]. The coupling of the leads and the central systems in the continuous representation conserves parity of the electron states across the tunneling barrier.

The full strength of the continuous approach emerges as it is applied to describe the transport of interacting electrons through 3D photon cavities in the transient time regime or the long time regime ending in a steady state of the system. This will be reported below (see Sections 5 and 6). The numerical calculations can sometimes be simplified by describing the leads and the central system on a discrete spatial lattice, where the geometric details of the central system are usually implemented by hard walls and Dirichlet boundary conditions. The spatial integral of the coupling tensor (37) are then reduced to a set of contact points between the leads and the central system [12,17].

### **3. Many-Body Effects in the Transient Regime**

In this section we review some results on the transient transport in interacting systems described by a lattice model [12,14]. For the sake of generality we extend the GME method by including as well the spin degree or freedom which was previously neglected. The lattice model matches naturally to the partitioning transport setting, facilitates the geometrical description of the central sample (e.g., a parallel quantum dot) and captures the dependence of the tunneling coefficients on the localization of the single-particle wavefunctions at the contact regions. A more realistic description is provided by the continuous model (see the previous section) which requires however a very careful tailoring of the confining potentials.

The results presented in this section are also meant to illustrate the usefulness of the GME approach in describing the transient regime in terms of the dynamical occupations of the interacting many-body configurations. Such a description cannot be recovered within the non-equilibrium Greens' function formalism.

Developing the GME method in the language of interacting many-body states was equally motivated by experimental works. Recording the charging of excited states of QDs in the Coulomb blockade regime constitutes the core of transient current spectroscopy and pump-and-probe techniques [54]. Also, transient currents through split-gate quantum point contacts (QPSs) and Ge quantum dots have been measured some time ago by Nasser et al. [55] and by Lai et al. [56]. Another relevant class of transport phenomena which can be modeled and understood within the GME method is the electron pumping through QDs with tunable-barriers (see e.g., the recent review [57]). In this context we investigated the transient response of a quantum dot submitted to a sequence of rectangular pulses applied at the contact to the input [58] and the turnstile protocol for single-molecule magnets [59].
