Generalized Master Equation Approach to Time-Dependent Many-Body Transport

Abstract

We recall theoretical studies on transient transport through interacting mesoscopic systems. It is shown that a generalized master equation (GME) written and solved in terms of many-body states provides the suitable formal framework to capture both the effects of the Coulomb interaction and electron–photon coupling due to a surrounding single-mode cavity. We outline the derivation of this equation within the Nakajima–Zwanzig formalism and point out technical problems related to its numerical implementation for more realistic systems which can neither be described by non-interacting two-level models nor by a steady-state Markov–Lindblad equation. We first solve the GME for a lattice model and discuss the dynamics of many-body states in a two-dimensional nanowire, the dynamical onset of the current-current correlations in electrostatically coupled parallel quantum dots and transient thermoelectric properties. Secondly, we rely on a continuous model to get the Rabi oscillations of the photocurrent through a double-dot etched in a nanowire and embedded in a quantum cavity. A many-body Markovian version of the GME for cavity-coupled systems is also presented.

1. Introduction

Few-level open systems stand as everyday ‘lab rats’ in corner stone experiments and future technologies in nanoelectronics [1] and quantum optics [2]. Generically, they are electronic systems with a discrete spectrum (e.g., artificial atoms [3], nanowires or superconducting qubits [4]) connected to particle reservoirs or embedded in bosonic baths. Depending on the nature of the environment (i.e., fermionic or bosonic) to which the open systems are coupled, their theoretical investigation started with two toy-models, namely the single-level Hamiltonian of quantum transport and the Jaynes–Cummings (JC) Hamiltonian of a two-level system (TLS).

Surprisingly or not, studying the sequential tunneling transport regime or the optical properties of quantum emitters eventually boils down to solve formally similar Markovian master equations (MEs) for the so called reduced density operator (RDO). The latter defines the non-unitary evolution of the small system in the presence of the infinite degrees of freedom of the reservoirs. Such MEs are derived by tracing out the reservoir’s degrees of freedom and are known from the early days of condensed matter and quantum optics (see the seminal works of Bloch [5], Wangsness [6] and Redfield [7]). The master equation cleverly bypasses the fact that the Liouville–von-Neumann (LvN) equation of the coupled systems (i.e., the open system and the reservoirs) is impossible to solve and takes advantage of the fact that all observables associated to the small and open system can be calculated as statistical averages w.r.t. the RDO.

Indeed, the RDO associated to the Jaynes–Cummings model has been a central object in quantum optics [8,9] (e.g., in the study of lasing and for the calculation of photon correlation functions). In this context the master equation (ME) approach goes as follows: (i) one studies an atomic few-level system whose eigenvalues and eigenfunctions are supposed to be known; (ii) the dissipation in the system (e.g., cavity losses or various non-radiative recombination processes) is included through the so called ‘jump’ operators; (iii) the occupation of atomic levels changes due to photon emission or absorption, but the particle number is conserved as the system is not coupled to particle reservoirs; (iv) under the Markov approximation the ME acquires a Lindblad form, usually solved in the steady-state regime.

The above scenario changes when one aims to derive a quantum master equation describing transport phenomena. (i) The Coulomb interaction effects on the spectrum and eigenstates of the system cannot be always neglected, especially for confined systems like quantum dots or nanowires; this requires a many-body derivation of the master equation; (ii) The tunneling between source/drain probes prevents the charge conservation in the central system and the main quantity of interest is the electronic current; (iii) Finally, the steady-state regime does not cover the whole physics and cannot even be guaranteed in general; moreover, the validity of the rotating-wave (RWA) and Markov approximations must be established more carefully [10,11]. In fact it turns out that when applied to transport processes the master equation must rather be solved in its non-Markovian version.

Such generalized master equations which take into account the memory effects have been mostly derived and implemented for time-dependent transport in non-interacting [12,13] and interacting [14] quantum dots, nanowires, and rings [15]. It turns out that the generalized master equation (GME) method is a valuable tool for modeling and monitoring the dynamics of specific many-body states as well as for investigating time-dependent propagation along a sample [16] or capturing charge sensing effects [17] and counting statistics in electrostatically parallel QDs [18]. In particular, Harbola et al. [19] showed that a Lindblad form of the quantum master equation is still recovered in the high bias limit and by assuming the RWA.

Since then, a lot of theoretical work has been done to improve and refine the quantum master equation formalism. A formally exact memory-kernel for the Anderson model was derived and calculated using real-time path integral Monte Carlo methods [20]. A hierarchical quantum master equation approach with good convergence at not too low temperature was put forward by Härtle et al. [21]. As for molecular transport calculations one can rely on the GME written in terms of the many-body states of the isolated molecule [22,23]. A recent review on non-markovian effects in open systems is also available [24].

As we shall see below the implementation of GME approach to many-level systems with specific geometries poses considerable technical difficulties. These are related to the many-body structure of the central interacting system, to the accurate description of the contact regions and, more importantly, to the evaluation of the non-Markovian kernels which become complicated objects once we go beyond non-interacting single-level models.

A second useful extension of the ME method emerged in the context of cavity quantum electrodynamics. Here the system under study is a hybrid one, as the electronic system is still coupled to source/drain reservoirs (i.e., leads) but also interacts with a quantum cavity mode, the latter being subjected to dissipation into leaky modes described by a bosonic bath. Such systems are currently used in state-of-the-art measurements in cavity quantum electrodynamics [25,26,27,28,29]. Again, the many-body nature of the problem is essential, as the electron-photon coupling leads to the formation of dressed states whose dynamics in the presence of both particle and dissipative bosonic reservoirs is far from being trivial. The relevant reduced density operator now acts in the many-body electron-photon Fock space and describes the dynamics of dressed-states. This fact brings new technical difficulties in the derivation [30,31] and implementation of ME [32,33]. Let us also mention here recent studies on ground state electroluminescence [34,35] and on cavity enhanced transport of charge [36].

In view of the abovementioned comments, the aim of this work is: (i) to briefly review the development of the generalized master equation approach to time-dependent many-body transport in the presence of both fermionic and bosonic environments and (ii) to illustrate in a unified framework how the method really works, from formal technicalities to numerical implementation. In Section 2 we shall therefore derive a non-Markovian master equation which describes the dynamics and the transport properties of rather general ‘hybrid’ system consisting in an electronic component $S_1$ which is connected to particle reservoirs (i.e., leads) and a second subsystem $S_2$. The latter, although not coupled to particle reservoirs, interacts with system $S_1$ or with some leaky modes described as bosonic baths. Then we specialize this master equation to several systems of interest. More precisely, in Section 3 we recall GME results on transient charging of excited states and Coulomb-coupled quantum dots. Section 4 deals with thermoelectric transport. Applications to transport in cavity quantum electrodynamics are collected in Section 5 and Section 6. We conclude in Section 7.

2. Formalism

2.1. Generalized Master Equation for Hybrid Systems

Non-Markovian master equations for open systems have been derived in many recent textbooks or review papers via projection methods (e.g., Nakajima–Zwanzig formalism or time-convolutionless approach [37]). Nonetheless it is still instructive to outline here some theoretical and computational difficulties one encounters when solving transport master equation for interacting many-level systems.

From the formal point of view the projection technique is quite general and the derivation of a master equation for the RDO does not depend on a specific model (i.e., on the geometry and spectrum of the central system or on the correlation functions of fermionic/bosonic reservoirs). In general, as long as one can write down a system-reservoir coupling Hamiltonian $H_{SR}$ a master equation can be derived.

For the sake of generality we shall consider a hybrid system S made of an electronic structure $S_1$ which is coupled to $n_r$ particle reservoirs characterized by chemical potentials and temperatures $\{{\mu }_l,T_l\}$, $l=1,2,\dots ,n_r$, and a second subsystem $S_2$ (i.e., a localized impurity, or an oscillator, or a single-mode quantum cavity). The subsystem $S_2$ can only be coupled to thermal or photonic baths which are described as a collection of oscillators with frequencies $\left\{{\omega }_k\right\}$. Let ${\mathcal{F}}_{S_1}$ and ${\mathcal{F}}_{S_2}$ be the Fock spaces associated to the two systems. Typically ${\mathcal{F}}_{S_1}$ is a set of interacting many-body configurations of the electronic system whereas ${\mathcal{F}}_{S_2}$ is made by harmonic oscillator Fock states.

The dynamics of the open system $S_1$ and of nearby ‘detector’ system $S_2$ are intertwined by a coupling V. Under a voltage bias or a temperature difference the system $S_1$ carries an electronic or a heat current which need to be calculated in the presence of the second subsystem. Conversely, the averaged observables of $S_2$ (e.g., mean photon number or the spin of a localized impurity) will also depend on the transport properties of $S_1$. The Hamiltonian of the hybrid structure is:

$$H_S=H_{S_1}+H_{S_2}+V.$$

(1)

In this work $H_{S_1}$ will describe various Coulomb-interacting structures: a single quantum dot, a 2D wire or parallel quantum dots. We shall denote by $|\nu \rangle$ and $E_{\nu }$ the many-body configurations and eigenvalues of $H_{S_1}$, that is one has $H_{S_1}|\nu \rangle =E_{\nu }|\nu \rangle$. $H_{S_1}$ can be equally expressed in terms of creation and annihilation operators $\{c_{n\sigma }^{†},c_{n\sigma }\}$ associated to a spin-dependent single-particle basis $\left\{{\psi }_{n\sigma }\right\}$ of a single-particle Hamiltonian $h_{S_1}^{\left(0\right)}$ (see the next sections for specific models), such that:

$$H_{S_1}=H_{S_1}^{\left(0\right)}+W,$$

(2)

where $H_{S_1}^{\left(0\right)}$ is the 2nd quantized form of $h_{S_1}^{\left(0\right)}$ and W is the Coulomb interaction. Similarly, the eigenstates and eigenvalues of the second subsystem $S_2$ will be denoted by $|j\rangle$ and $e_j$ such that $H_{S_2}|j\rangle =e_j|j\rangle$. As for the coupling V one can mention at least three examples: The exchange interaction between a quantum dot and a localized magnetic impurity with total spin S, the electron–photon coupling in a quantum-dot cavity and the electron–vibron coupling in nanoelectromechanical systems [38,39].

The total Hamiltonian of the system coupled to particle and/or bosonic reservoirs R reads as:

$$H\left(t\right)=H_S+H_R+H_{SR}\left(t\right):=H_0+H_{SR}\left(t\right),$$

(3)

where the system-reservoir coupling $H_{SR}$ collects the coupling to the leads ($H_T$) and the coupling of a bosonic mode to a thermal or leaky bosonic environments ($H_E$):

$$H_{SR}\left(t\right)=H_T\left(t\right)+H_E.$$

(4)

Note that the interaction with the bosonic environment $H_E$ does not depend on time. The lead-sample tunneling term $H_T$ carries a time-dependence that will be explained below. The Hamiltonian of the reservoirs,

$$H_R=H_{\mathrm{leads}}+H_{\mathrm{bath}}$$

(5)

describes at least two semiinfinite leads (left-L and right-R) but could also contain a bosonic or a thermal bath.

This general scheme allows one to recover several relevant settings. If $S_1$ describes an optically active structure and $S_2$ defines a photonic mode then V could become either the Rabi or the Jaynes–Cummings electron–photon coupling. The absence of the particle reservoirs simplifies $H_S$ to well known models in quantum optics, while by adding them one can study photon-assisted transport effects (e.g., Rabi oscillation of the photocurrents or electroluminescence). Also, by removing $S_2$, V and the bosonic dissipation one finds the usual transport setting for a Coulomb interacting purely electronic structure.

Let ${\epsilon }^l\left(q\right)$ and ${\psi }_{q\sigma }^l$ be the single particle energies and wave functions of the l-th lead. For simplicity we assume that the states on the leads are spin-degenerate so their energy levels do not depend on the spin index. Using the creation/annihilation operators $c_{ql\sigma }^{†}$/$c_{ql\sigma }$ associated to the single particle states, we can write:

$$H_{\mathrm{leads}}=\sum _l{H}_l=\int dq\sum _{\sigma }{\epsilon }^l\left(q\right)c_{ql\sigma }^{†}{c}_{ql\sigma }.$$

(6)

As for the bosonic bath, it is described by a collection of harmonic oscillators with frequencies ${\omega }_k$ and by corresponding creation/annihilation operators $b_k^{†}/b_k$:

$$H_{\mathrm{bath}}=\sum _k\hslash {\omega }_k{b}_k^{†}{b}_k.$$

(7)

The tunneling Hamiltonian has the usual form:

$$H_T\left(t\right)=\sum _l\sum _{n\sigma }\int dq{\chi }_l\left(t\right)(T_{qn}^l{c}_{ql\sigma }^{†}{c}_{n\sigma }+h.c),$$

(8)

where we considered without loss of generality that the tunneling processes are spin conserving. For the simplicity of writing the spin degree of freedom $\sigma$ will be henceforth tacitly merged with the single-particle index n and restored if needed.

The time-dependent switching functions ${\chi }_l\left(t\right)$ control the time-dependence of the contacts between the leads and the sample; these functions mimic the presence of a time dependent potential barrier. We emphasize that in most studies based on ME method the coupling to the leads is suddenly switched at some initial instant $t_0$ such that for each lead ${\chi }_l\left(t\right)=\theta (t-t_0)$ where $\theta \left(x\right)$ is the Heaviside step function. This choice is very convenient if one imposes the Markov approximation in view of a time-local Master equation. Here we allow for more general switching functions: (i) a smooth coupling to the leads or (ii) time-dependent signals applied at the contacts to the leads. In particular, if the potential barriers oscillate out of phase the system operates like a turnstile pump under a finite constant bias.

The coupling $T_{qn}^l$ describes the tunneling strength between a state with momentum q of the lead l and the state n of the isolated sample with wavefunctions ${\psi }_n$. In the next sections we shall show that these matrix elements have to be calculated for each specific model by taking into account the geometry of the system and of the leads.

The associated density operator $\mathcal{W}$ of the open system obeys the Liouville–von Neumann equation:

$$i\hslash \frac{\partial \mathcal{W}\left(t\right)}{\partial t}=\mathcal{L}\left(t\right)\mathcal{W}\left(t\right),\mathcal{W}\left(t_0\right)={\rho }_S\left(t_0\right)\otimes {\rho }_R,$$

(9)

where:

$$\mathcal{L}\left(t\right)={\mathcal{L}}_0+{\mathcal{L}}_{SR},{\mathcal{L}}_0·=[H_0,·].$$

(10)

We also introduce the notations:

$${\mathcal{L}}_S·=[H_S,·],{\mathcal{L}}_{SR}·=[H_{SR},·].$$

(11)

The Nakajima–Zwanzig projection formalism leads to an equation of motion for the reduced density operator $\rho \left(t\right)={\mathrm{Tr}}_R\left\{\mathcal{W}\right\}$. The initial state ${\mathcal{W}}_0:=\mathcal{W}\left(t_0\right)$ factorizes as:

$${\mathcal{W}}_0={\rho }_0\otimes {\rho }_{\mathrm{leads}}\otimes {\rho }_{\mathrm{bath}}:={\rho }_0\otimes {\rho }_R,$$

(12)

where the equilibrium density operator of the leads reads:

$${\rho }_{\mathrm{leads}}=\prod _l\frac{e^{-{\beta }_l(H_l-{\mu }_l{N}_l)}}{{\mathrm{Tr}}_l\left\{e^{-{\beta }_l(H_l-{\mu }_l{N}_l)}\right\}},$$

(13)

and ${\beta }_l=1/k_B{T}_l$, ${\mu }_l$ and $N_l$ denote the inverse temperature, chemical potential and the occupation number operator of the lead l. Similarly,

$${\rho }_{\mathrm{bath}}=\prod _k{e}^{-\hslash {\omega }_k{b}_k^{†}{b}_k/k_{B}T}(1-e^{-\hslash {\omega }_k/k_{B}T}).$$

(14)

Finally, ${\rho }_0$ is simply a projection on one of the states of the hybrid system, and as such its calculation must take into account the effect of the hybrid coupling V (see the discussion in Section 2.2). We now define two projections:

$$P·={\rho }_R{\mathrm{Tr}}_R\{·\},Q=1-P.$$

(15)

It is straightforward to check the following properties:

$$\begin{array}{ccc}\hfill P{\mathcal{L}}_S& =& {\mathcal{L}}_{S}P,P{\mathcal{L}}_{SR}\left(t\right)P=0.\hfill \end{array}$$

(16)

The Liouville Equation (9) splits then into two equations:

$$\begin{array}{ccc}\hfill i\hslash P\frac{\partial W\left(t\right)}{\partial t}& =& P\mathcal{L}\left(t\right)PW\left(t\right)+P\mathcal{L}\left(t\right)QW\left(t\right)\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill i\hslash Q\frac{\partial W\left(t\right)}{\partial t}& =& Q\mathcal{L}\left(t\right)QW\left(t\right)+Q\mathcal{L}\left(t\right)PW\left(t\right),\hfill \end{array}$$

(18)

and the second equation can be solved by iterations (T being the time-ordering operator):

$$QW\left(t\right)=\frac{1}{i\hslash }{\int }_{t_0}^{t}dsTexp\left\{-\frac{i}{\hslash }{\int }_s^{t}d{s}^{\prime }Q\mathcal{L}\left(s^{\prime }\right)\right\}Q\mathcal{L}\left(s\right)PW\left(s\right).$$

(19)

Inserting Equation (19) in Equation (17) and using the properties of P we get the Nakajima–Zwanzig equation:

$$\begin{array}{ccc}\hfill i\hslash P\frac{\partial W\left(t\right)}{\partial t}& =& P{\mathcal{L}}_{S}W\left(t\right)\hfill \\ \hfill & +& \frac{1}{i\hslash }P{\mathcal{L}}_{SR}\left(t\right)Q{\int }_{t_0}^{t}dsTexp\left\{-\frac{i}{\hslash }{\int }_s^{t}d{s}^{\prime }Q\mathcal{L}\left(s^{\prime }\right)Q\right\}Q{\mathcal{L}}_{SR}\left(s\right)PW\left(s\right).\hfill \end{array}$$

(20)

In order to have an explicit perturbative expansion in powers of $H_{SR}\left(t\right)$ one has to factorize the time-ordered exponential as follows:

$$Texp\left\{-\frac{i}{\hslash }{\int }_s^{t}d{s}^{\prime }Q\mathcal{L}\left(s^{\prime }\right)Q\right\}=exp\left\{Q{\mathcal{L}}_{0}Q\right\}(1+\mathcal{R}),$$

(21)

where the remainder $\mathcal{R}$ contains infinitely deep commutators with inconveniently embedded projection operators. Usually one considers a truncated version of the Nakajima–Zwanzig equation up to the second order contribution w.r.t. the system-reservoir $H_{SR}$:

$$i\hslash \dot{\rho }\left(t\right)={\mathcal{L}}_S\rho \left(t\right)+\frac{1}{i\hslash }{\mathrm{Tr}}_R\left\{{\mathcal{L}}_{SR}{\int }_{t_0}^{t}ds{e}^{-i(t-s){\mathcal{L}}_0}{\mathcal{L}}_{SR}\left(s\right){\rho }_R\rho \left(s\right)\right\}.$$

(22)

Now, by taking into account that for any operator A acting on the Fock space of the hybrid system $e^{-i{t\mathcal{L}}_0}A=e^{-it{H}_0}A{e}^{it{H}_0}$ and denoting by $U_0(t,s)=e^{-i(t-s)H_0}$ the unitary evolution of the disconnected systems we arrive at the well known form of the GME:

$$\begin{array}{ccc}i\hslash \dot{\rho }\left(t\right)\hfill & =\hfill & [H_S,\rho \left(t\right)]-\frac{i}{\hslash }{U}_0^{†}(t,t_0){\mathrm{Tr}}_R\left\{{\int }_{t_0}^{t}ds\left[{\tilde{H}}_{SR}\left(t\right),\left[{\tilde{H}}_{SR}\left(s\right),\tilde{\rho }\left(s\right){\rho }_R\right]\right]\right\}{U}_0(t,t_0)\hfill \\ \hfill & =\hfill & [H_S,\rho \left(t\right)]-\frac{i}{\hslash }{U}_S^{†}(t,t_0){\mathrm{Tr}}_R\left\{{\int }_{t_0}^{t}ds\left[{\tilde{H}}_{SR}\left(t\right),\left[{\tilde{H}}_{SR}\left(s\right),\tilde{\rho }\left(s\right){\rho }_R\right]\right]\right\}{U}_S(t,t_0),\hfill \end{array}$$

where in order to get to the last line we removed the evolution operators of the environment from both sides of the trace. At the next step one observes that when performing the trace over the reservoirs and environment degrees of freedom the mixed terms in the double commutator vanish because each of the coupling terms $H_T$ and $H_E$ carries only one creation or annihilation operator for the corresponding reservoir such that:

$${\mathrm{Tr}}_R\left\{{\tilde{c}}_{ql}^{†}\left(t\right){\tilde{b}}_k\left(s\right){\rho }_R\right\}={\mathrm{Tr}}_{\mathrm{leads}}\left\{{\tilde{c}}_{ql}^{†}\left(t\right){\rho }_{\mathrm{leads}}\right\}·{\mathrm{Tr}}_{\mathrm{bath}}\left\{{\tilde{b}}_k\left(s\right){\rho }_{\mathrm{bath}}\right\}=0.$$

(23)

Moreover, the time evolution of each term can be simplified due to the commutation relations $[H_{\mathrm{bath}},H_T]=[H_{\mathrm{leads}},H_E]=0$:

$$\begin{array}{ccc}\hfill {\tilde{H}}_T\left(t\right)& =& e^{\frac{i}{\hslash }t{H}_S}{e}^{\frac{i}{\hslash }t{H}_{\mathrm{leads}}}{H}_T{e}^{-\frac{i}{\hslash }t{H}_S}{e}^{-\frac{i}{\hslash }t{H}_{\mathrm{leads}}},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\tilde{H}}_E\left(t\right)& =& e^{\frac{i}{\hslash }t{H}_S}{e}^{\frac{i}{\hslash }t{H}_{\mathrm{bath}}}{H}_E{e}^{-\frac{i}{\hslash }t{H}_S}{e}^{-\frac{i}{\hslash }t{H}_{\mathrm{bath}}}.\hfill \end{array}$$

(25)

The GME then reads as:

$$\begin{array}{ccc}\hfill \dot{\rho }\left(t\right)& =& -\frac{i}{\hslash }[H_S,\rho \left(t\right)]-\frac{1}{{\hslash }^2}{U}_S^{†}(t,t_0){\mathrm{Tr}}_{\mathrm{leads}}\left\{{\int }_{t_0}^{t}ds\left[{\tilde{H}}_T\left(t\right),\left[{\tilde{H}}_T\left(s\right),\tilde{\rho }\left(s\right){\rho }_{\mathrm{leads}}\right]\right]\right\}{U}_S(t,t_0)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill & -& \frac{1}{{\hslash }^2}{U}_S^{†}(t,t_0){\mathrm{Tr}}_{\mathrm{bath}}\left\{{\int }_{t_0}^{t}ds\left[{\tilde{H}}_E\left(t\right),\left[{\tilde{H}}_E\left(s\right),\tilde{\rho }\left(s\right){\rho }_{\mathrm{bath}}\right]\right]\right\}{U}_S(t,t_0)\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill & :=& -\frac{i}{\hslash }[H_S,\rho \left(t\right)]-{\mathcal{D}}_{\mathrm{leads}}[\rho ,t]-{\mathcal{D}}_{\mathrm{bath}}[\rho ,t].\hfill \end{array}$$

(27)

Equation (27) is the generalized master equation for our hybrid system. It provides the reduced density operator $\rho$ in the presence of particle and bosonic reservoirs and also takes into account the memory effects and the non-trivial role of time-dependent signals applied at the contact regions through the switching functions ${\chi }_l$. The third term in Equation (27) is needed only if $H_{S_2}$ describes a quantized optical or mechanical oscillation mode. In our work on open QD-cavity systems we always assume that the coupling between the cavity photons and other leaky modes is much smaller that the electron-photon coupling $g_{EM}$ (see Section 5). On the other hand, for our calculations $g_{EM}\ll \hslash \omega$, $\omega$ being the frequency of the cavity mode. Then the RWA holds and ${\mathcal{D}}_{\mathrm{bath}}[\rho ,t]$ can be cast in a Lindblad form. Let us stress that in the ultrastrong coupling regime on typically has $g_{EM}/\hslash \omega >0.2$ and the derivation of the dissipative term is more complicated and involves the dressed states of the QD-cavity system [30]. In order to describe dissipative effects in the ultrastrong coupling regime beyond the RWA one needs more elaborate techniques [40,41].

For further calculations one has to solve the GME as a system of coupled integro-differential equations for the matrix elements of the RDO with respect to a suitable basis in the Fock space ${\mathcal{F}}_S={\mathcal{F}}_{S_1}\otimes {\mathcal{F}}_{S_2}$. We discuss this issue in the next subsection.

2.2. ‘Hybrid’ States and Diagonalization Procedure

The starting point in solving the GME is to write down the matrix elements of the system-environment operators $H_T$ and $H_E$ w.r.t. the ‘disjointed’ basis formed by the eigenstates of $H_{S_1}$ and $H_{S_2}$, that is $|\nu ,j\rangle :=|\nu \rangle \otimes |j\rangle$. However this strategy does not help much when evaluating the time evolution (see Equations (24) and (25)) as $H_S$ is not diagonal w.r.t. to $|\nu ,j\rangle$ such that one cannot easily write down the matrix elements of the unitary evolution $U_S(t,t_0)$. In fact we are forced to solve the GME by using the eigenstates $|{\phi }_p)$ and eigenvalues ${\mathcal{E}}_p$ of the Hamiltonian $H_S$. The former are written as:

$$|{\phi }_p)=\sum _{\nu ,j}{\mathcal{V}}_{\nu j}^{\left(p\right)}|\nu ,j\rangle .$$

(28)

Here the round bracket notation $|{\phi }_p)$ is meant to underline that the state ${\phi }_p$ describes the interacting system S, in the sense that both Coulomb interactions and the coupling to the bosonic modes were taken into account when diagonalizing $H_S$. This notation also prevents any confusion if the ‘free’ states $|\nu ,j\rangle$ were also labelled by a single index $p^{\prime }$. In that case the above equation is conveniently rewritten as $|{\phi }_p)={\sum }_{p^{\prime }}{\mathcal{V}}_{p^{\prime }}^{\left(p\right)}|p^{\prime }\rangle$. Note that p is usually a multiindex carrying information on relevant quantum numbers. In most cases of interest the coupling V between the two systems leads to a strong mixing of the unperturbed basis elements $|\nu ,j\rangle$ and is not necessarily small. Therefore we shall not follow a perturbative approach but rather calculate ${\mathcal{E}}_p$ and the weights ${\mathcal{V}}_{\nu j}^{\left(p\right)}$ by numerically diagonalizing $H_S$ on a relevant subspace of ‘disjointed’ states.

Prior to any model specific calculations or numerical implementations it is useful to comment a bit on the two dissipative contributions in Equation (27). It is clear that the evolution operator $U_S$ describes the joint systems $S_1$ and $S_2$ and therefore the hybrid interaction cannot be simply neglected neither in ${\mathcal{D}}_{\mathrm{leads}}$ nor in ${\mathcal{D}}_{\mathrm{bath}}$; in fact one can easily check that V does not commute with $H_E$ or $H_T$. Moreover, as has been clearly pointed out by Beaudoin et al. [30], by disregarding the qubit-resonator interaction when calculating ${\mathcal{D}}_{\mathrm{bath}}$ one ends up with unphysical results. In what concerns the role of V in the leads’ contribution, a recent work emphasized that for QD-cavity systems the corresponding master equation must be derived in the basis of dressed-states [31].

The diagonalization of $H_S$ poses serious technical problems because both spaces ${\mathcal{F}}_{S_1}$ and ${\mathcal{F}}_{S_2}$ are in principle infinite dimensional. Besides that, the Coulomb interaction in $H_{S_1}$ prevents one to derive the interacting many-body configurations $\{|\nu \rangle \}$ analytically. We now propose a step-by-step diagonalization procedure leading to a relevant set of interacting states of the full Hamiltonian. The procedure requires several ‘intermediate’ diagonalization operations:

($D_1$) Analytical or numerical calculation of the single-particle states of the Hamiltonian ${\widehat{h}}_{S_1}^{\left(0\right)}$ which describes the non-interacting electronic system $S_1$. As we shall see in the next sections, this step may not be trivial if the geometry of the sample is taken seriously into account. Let us select a subset of $N_{\mathrm{ses}}$ single-particle states $\{{\psi }_1,{\psi }_2,\dots ,{\psi }_{N_{\mathrm{ses}}}\}$ (if needed this set of states includes the spin degree of freedom). Typically we choose the lowest energy single-particle states but in some cases [12] it is more appropriate to select the subset of states which effectively contribute to the transport (i.e., states located within the bias window).

($D_2$) The construction of a second set of $N_{\mathrm{mes}}$ non-interacting many-body configurations (NMBS) ${\{|\lambda \rangle \}}_{\lambda =1,\dots ,N_{\mathrm{mes}}}$ from the $N_{\mathrm{ses}}$ single-particle states introduced above. Note that for computational reasons we have to keep $N_{\mathrm{mes}}<2^{N_{\mathrm{ses}}}$ for larger $N_{\mathrm{ses}}$. Then, if $i_n^{\lambda }$ is the occupation number of the single-particle state ${\psi }_n$, a non-interacting many-body configuration $|\lambda \rangle$ reads as:

$$|\lambda \rangle =|i_1^{\lambda },i_1^{\lambda },\dots ,i_{N_{\mathrm{ses}}}^{\lambda }\rangle .$$

(29)

($D_3$) Diagonalization of the Coulomb-interacting electronic Hamiltonian $H_{S_1}=H_{S_1}^{\left(0\right)}+W$ on the subspace of non-interacting many-body states from ${\mathcal{F}}_{S_1}$. As a result one gets $N_{\mathrm{mes}}$ interacting many-body states (IMBS) $|\nu \rangle$ and the associated energy levels $E_{\nu }$ introduced in Section 2.1. We also introduce the ‘free’ energies of $H_{S_1}+H_{S_2}$, that is ${\mathcal{E}}_{\nu ,j}^{\left(0\right)}=E_{\nu }+e_j$. Note that in view of diagonalization the interaction V between the two subsystems must be also written w.r.t. the ‘free’ states $\{|\nu ,j\rangle \}$. If the second subsystem is not needed then the GME must be solved w.r.t. the set $\{|\nu \rangle \}$ and the diagonalization procedure stops here. It is worth pointing out here that even in the absence of bosonic fields and electron-photon coupling, the master equation for Coulomb interacting systems cannot be written in terms of single particle states. In spite of the fact that the unitary evolution $U_S$ is diagonal w.r.t. the many-body basis $\{|\nu \rangle \}$, the matrix elements of the fermionic operators in the interaction picture $\langle \nu |{\tilde{c}}_{n\sigma }\left(t\right)|{\nu }^{\prime }\rangle$ depend on energy differences $E_{\nu }-E_{{\nu }^{\prime }}$ which, due to the Coulomb effects, cannot be reduced to the single-particle energy ${\epsilon }_n$.

($D_4$) Diagonalization of the fully interacting Hamiltonian $H_S$ on a subspace of $\mathcal{F}$ made by the lowest energy $N_{\mathrm{mesT}}$ interacting MBS of $H_{S_1}$ and $j_{\max }$ eigenstates of $H_{S_2}$. Remark that after the 1st truncation w.r.t. NMBSs we perform here a 2nd double truncation both w.r.t. IMBS (as $N_{\mathrm{mesT}}<N_{\mathrm{mes}}$) and w.r.t. the states in ${\mathcal{F}}_{S_2}$.

Once this procedure is performed, one can express the system-environment couplings $H_T$ and $H_B$ in the fully interacting basis and use the eigenvalues ${\mathcal{E}}_p$ to replace the unitary evolution $U_S$ by the corresponding diagonal matrix $e^{-it{\mathcal{E}}_p}{\delta }_{p{p}^{\prime }}$. Finally, the GME is to be solved w.r.t. the fully interacting basis (see Section 2.3).

Now let us enumerate and explain the advantages of this stepwise procedure when compared to a single and direct diagonalization of $H_S$.

(a) Numerical efficiency and accuracy. Both diagonalization methods (stepwise and direct) require a truncation of both bases and are not free of numerical errors which in principle should diminish as the size of the bases increase. It is clear that in the stepwise procedure the $N_{\mathrm{mes}T}$ interacting MBSs are derived from a larger set of non-interacting states ${\{|\lambda \rangle \}}_{\lambda =1,\dots ,N_{\mathrm{mes}}}$. Then the calculated $N_p:=N_{\mathrm{mes}T}\times j_{\max }$ fully interacting states are more accurate than the ones obtained by diagonalizing once a $N_p\times N_p$ matrix. On the other hand, enlarging the full space to $N_{\mathrm{mes}}\times j_{\max }$ elements could be challenging in terms of CPU times. Convergence calculations relevant to circuit quantum electrodynamics have been presented in [42]. In particular it was shown that the inclusion of the (usually neglected) diamagnetic term in the electron-photon coupling improves the convergence of the diagonalization procedure.

The size of various effective bases used in the numerical calculation is decided both by physical considerations and convergence tests. Typically, out of the $N_{\mathrm{ses}}$ single-electron states we construct the set of non-interacting MBSs containing up to $N_e$ electrons, the size of this set being, of course, $\left(\genfrac{}{}{0pt}{}{N_{\mathrm{ses}}}{N_e}\right)$. The accuracy of numerical diagonalization which leads to the interacting many-body configurations with up to $N_e$ electrons is essentially assessed by comparing the spectra associated obtained for different $N_{\mathrm{ses}}$. In particular, if we discretize our open system in a small number of lattice sites we can use all single-electron states as a basis, and we can calculate all many-body electron states (like in the discrete case presented in [12]. Obviously, this is no longer possible for a more complex geometry, and then we need to evaluate the convergence of the results when the basis is truncated.

For the continuous model an extensive discussion on the convergence of the numerical diagonalization w.r.t. the various truncated bases is given in a previous publication [42]. Let us stresonly extensive tests can be performed to resolve this issue. s here that once the geometry of the system and the spatially-dependent interactions are accounted for there is no simple way to count ahead how many states one needs to get stable transport simulations.

(b) Physical interpretation. It is obvious that the Coulomb interaction W mixes only the non-interacting many-body configurations $|\lambda \rangle$ while the hybrid coupling V mixes both $\lambda$ and j states. For this reason the weights of a non-interacting state $|\lambda ,j\rangle$ in a fully interacting state $|{\phi }_p)$ (as provided by a single diagonalization) cannot be easily associated with one of the interacting terms. In view of physical discussion it is more intuitive and natural to analyze the dynamics of the Coulomb-interacting system $S_1$ in the presence of the second subsystem $S_2$. One such example is a self-assembled quantum dot embedded in a single-mode quantum cavity [31]. In this system the optical transitions couple electron-hole pairs which are genuine interacting many-body states. A second example is a double quantum dot patterned in a 2D quantum wire which is itself placed in a cavity. There the interdot Coulomb interaction affects the optical transitions as well.

On the other hand, the above procedure will not be appropriate if one is interested in including a time-dependent driving term in $H_S$. This would be the case for a pumping potential or for a modulating optical signal.

Finally we shall comment on the typical sizes of many-body Fock spaces used to model steady-state or transient transport in various systems. One of the simplest yet promising system for solid-state quantum computation is a double quantum dot accomodating at most two electrons on each dot such that the relevant Fock space already comprises 256 interacting many-body states (counting the spin). In this case transport simulations can be obtained even without truncating the basis, especially if the spectral gaps allows one to disregard the contribution of higher energy configurations. However, when studying transport on edge states due to a strong perpendicular magnetic field in 2D systems (e.g., graphene or phosforene) one is forced to consider the low energy bulk-states as well. In our previous numerical studies we find that one has to take into account at least 10 single-particle states; obviously, performing time-dependent simulations for $2^{10}$ MBSs is quite inconvenient so a truncation is needed. More importantly, a realistic description of complex systems like rings or double dots etched in a 2DEG cannot be obtained with only few single-particle states. Note that the optical selection rules and matrix elements of the electron–photon interaction depend on the these states as well.

In the calculation of one, rather deep, QD embedded in a short quantum wire we are using 52 single-electron states, asking for 52 one-electron states, 1326 two-electron states and 560 three electron states. Of these we take the lowest in energy 512 and tensor multiply by 17 photon states to obtain a basis of 8704 MBS to calculate the dressed MBS. Then for the transport, we select the lowest in energy 128 dressed states and construct the 16,384 dimensional Liouville space. All this choice is taylored for a rather narrow section of a parameter space, if we consider the wire length, the confinement energy and the shape of the QD and the range of the magnetic field.

Markovian or non-markovian master equation method have been also developed for transport simulations in molecular junctions; here a truncation is required w.r.t. to the basis states describing the molecular vibrations. In particular, Schinabeck et al. [43] proposed a hierarchical polaron master equation which was successfully implemented numerically for two molecular orbitals and several tens of vibrational states.

2.3. Numerical Implementation and Observables

The last step before numerical implementation requires the calculation of the system-environment couplings $H_T$ and $H_E$ w.r.t. the full basis $|{\phi }_p)$. Clearly, to this end we shall use the unitary transformations $|\lambda \rangle \leftrightarrow |\nu \rangle$ and $|\nu ,j\rangle \leftrightarrow |{\phi }_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\left(q\right)=\sum _{p,p^{\prime }}{\mathcal{T}}_{p{p}^{\prime }}^l\left(q\right)|{\phi }_p){({\phi }_{p^{\prime }}|,{(\mathcal{T}}_l\left(q\right))}_{p{p}^{\prime }}=\sum _n{T}_{qn}^l({\phi }_p|c_n^{†}|{\phi }_{p^{\prime }}).$$

(30)

Then the dissipation operator associated to the particle reservoirs reads:

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

(31)

with the following notation:

$$\begin{array}{c}{\Omega }_{ql}\left(t\right)=U_{\mathrm{S}}(t,t_0){\int }_{t_0}^{t}ds{\chi }_l\left(s\right){\Pi }_{ql}\left(s\right)e^{i((s-t)/\hslash ){\epsilon }_l\left(q\right)}{U}_{\mathrm{S}}^{†}(t,t_0),\hfill \end{array}$$

(32)

$$\begin{array}{c}{\Pi }_{ql}\left(s\right)=U_{\mathrm{S}}^{†}(s,t_0)\left({\mathcal{T}}_l^{†}\rho \left(s\right)(1-f_l)-\rho \left(s\right){\mathcal{T}}_l^{†}{f}_l\right)U_{\mathrm{S}}(s,t_0),\hfill \end{array}$$

(33)

and where for simplicity we omit to write the energy dependence of the Fermi function $f_l$. Similarly, the bosonic operators have to written down w.r.t. the full basis which then leads to the calculation of ${\mathcal{D}}_{\mathrm{bath}}$. 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 $\rho \left(t_n\right)$, 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 ${\Omega }_{ql}\left(t_{n+1}\right)$ obeys a recursive formula generated by the incremental integration between $t_n$ and $t_{n+1}$, that is:

$${\Omega }_{ql}\left(t_{n+1}\right)=U_{\mathrm{S}}(t_{n+1},t_n){\Omega }_{ql}\left(t_n\right)U_{\mathrm{S}}^{†}(t_{n+1},t_n)+{\mathcal{A}}_{ql}(t_{n+1},t_n;\rho \left(t_{n+1}\right),\rho \left(t_n\right)),$$

(34)

where the second term of the right-hand side depends on the yet unknown $\rho \left(t_{n+1}\right)$. For any pair of time steps $\{t_n,t_{n+1}\}$ we initially approximate $\rho \left(t_{n+1}\right)$ in ${\mathcal{A}}_{ql}$ by the already calculated $\rho \left(t_n\right)$, and perform iterations to recalculate $\rho \left(t_{n+1}\right)$ via the GME, each time updating $\rho \left(t_{n+1}\right)$ in ${\mathcal{A}}_{ql}$, until a convergence test for $\rho \left(t_{n+1}\right)$ is fulfilled. At any step of the iteration we also calculate and include ${\mathcal{D}}_{\mathrm{bath}}[\rho ,t_m]$ 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 $\rho \left(t_m\right)$, i.e., the populations of fully interacting states $|{\phi }_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 $\mathcal{O}$ of the hybrid system:

$$\langle \mathcal{O}\rangle ={\mathrm{Tr}}_{\mathcal{F}}\left\{\rho \left(t\right)\mathcal{O}\right\}.$$

(35)

Useful examples are averages of the photon number operator ${\mathcal{N}}^{ph}=a^{†}a$ and of the charge operator $\mathcal{Q}={\sum }_n{c}_n^{†}{c}_n$. Also, the average currents in a two-lead geometry (i.e., $l=L,R$ can be identified from the continuity equation:

$$\langle \dot{\mathcal{Q}}\rangle ={\mathrm{Tr}}_{\mathcal{F}}\left\{\mathcal{Q}\dot{\rho }\left(t\right)\right\}=J_L\left(t\right)-J_R\left(t\right).$$

(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 $L_x$ 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 }_S^l\times {\Omega }_l}d\mathbf{r}d{\mathbf{r}}^{\prime }{\left({\Psi }_q^l\left({\mathbf{r}}^{\prime }\right)\right)}^{*}{\Psi }_n^S\left(\mathbf{r}\right)g_{qn}^l(\mathbf{r},{\mathbf{r}}^{\prime })+h.c.,$$

(37)

for states with wavefunction ${\Psi }_q^l$ in lead l, and ${\Psi }_n^S$ in the central system. The domains for the integration of the wavefunctions in the leads are chosen to be

$$\begin{array}{ccc}{\Omega }_{\mathrm{L}}& =& \left\{(x,y)|\left[-\frac{L_x}{2}-2a_w,-\frac{L_x}{2}\right]\times [-3a_w,+3a_w]\right\},\\ {\Omega }_{\mathrm{R}}& =& \left\{(x,y)|\left[+\frac{L_x}{2},+\frac{L_x}{2}+2a_w\right]\times [-3a_w,+3a_w]\right\},\end{array}$$

(38)

and for the system as

$$\begin{array}{ccc}{\Omega }_{\mathrm{S}}^{\mathrm{L}}& =& \left\{(x,y)|\left[-\frac{L_x}{2},-\frac{L_x}{2}+2a_w\right]\times [-3a_w,+3a_w]\right\},\\ {\Omega }_{\mathrm{S}}^{\mathrm{R}}& =& \left\{(x,y)|\left[+\frac{L_x}{2}-2a_w,+\frac{L_x}{2}\right]\times [-3a_w,+3a_w]\right\}.\end{array}$$

(39)

The function

$$g_{qn}^l(\mathbf{r},{\mathbf{r}}^{\prime })=g_0^{l}exp\left[-{\delta }_1^l{(x-x^{\prime })}^2-{\delta }_2^l{(y-y^{\prime })}^2\right]exp\left(\frac{-|E_n-{ϵ}^l\left(q\right)|}{{\Delta }_E^l}\right).$$

(40)

with $\mathbf{r}\in {\Omega }_{\mathrm{S}}^l$ and ${\mathbf{r}}^{\prime }\in {\Omega }_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 ${\delta }_1^l$ and ${\delta }_1^l$ define the spatial range of the coupling within the domains ${\Omega }_{\mathrm{S}}^l\times {\Omega }_l$ [44].

The short quantum wire is considered to have hard wall confinement in the transport direction, the x-direction, at $x=\pm L_x/2$, and parabolic confinement in the y-direction with characteristic energy $\hslash \Omega$. Possibly, the leads and the central system are considered to be placed in a perpendicular homogeneous external magnetic field $\mathbf{B}=B\widehat{\mathbf{z}}$. Together they lead to a natural length scale, the effective magnetic length $a_w$ with $a_w^2{\Omega }_w=\hslash /m$, with ${\Omega }_w^2={[{\left({\Omega }_0\right)}^2+{\left({\omega }_s\right)}^2]}^{1/2}$, and the cyclotron frequency ${\omega }_c=(eB/m)$. For GaAs with effective mass $m=0.067m_e$ relative dielectric constant ${\kappa }_c=12.3$ and confinement energy $\hslash {\Omega }_0=2.0$ meV, $a_w=23.8$ nm. The magnetic field $B=0.1$ T.

The energy spectrum of the quasi 1D semi-infinite lead l is represented by ${ϵ}^l\left(q\right)$, with q standing for the momentum quantum number of a standing wave state, and the subband index $n^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 Section 5 and Section 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].

3.1. Transient Charging of Excited States

We consider a two-dimensional system of length $L_x$ and width $L_y$ described by a lattice with $N_x$ sites on the x axis and $N_y$ sites on the y axis. The total number of sites is denoted by $N_{xy}=N_x{N}_y$. By setting the two lattice constants $a_x$ and $a_y$ one has $L_x=a_x{N}_x$ and $L_y=a_y{N}_y$. Once we know the single-particle eigenstates of the electronic subsytem $S_1$ we can write down its Hamiltonian $H_{S1}:=H_{S_1}^{\left(0\right)}+W$ in a second quantized form w.r.t. this basis, that is:

$$\begin{array}{ccc}{H}_{S_1}^{\left(0\right)}\hfill & =\hfill & {\sum }_{n,n^{\prime }=1}^{N_{xy}}\langle {\psi }_n|{\widehat{h}}_{S_1}^{\left(0\right)}|{\psi }_{n^{\prime }}\rangle c_n^{†}{c}_{n^{\prime }}={\sum }_n{ϵ}_n{c}_n^{†}{c}_n\hfill \\ W\hfill & =\hfill & \frac{1}{2}{\sum }_{n,m,n^{\prime },m^{\prime }}{V}_{nm{n}^{\prime }{m}^{\prime }}{c}_n^{†}{c}_m^{†}{c}_{m^{\prime }}{c}_{n^{\prime }}\hfill \end{array}$$

(41)

where the Coulomb matrix elements are given by ($\mathbf{r},{\mathbf{r}}^{\prime }$ are sites of the 2D lattice):

$$V_{nm{n}^{\prime }{m}^{\prime }}=\sum _{\mathbf{r},{\mathbf{r}}^{\prime }}{\psi }_n^{*}\left(\mathbf{r}\right){\psi }_m^{*}\left({\mathbf{r}}^{\prime }\right)V_C(\mathbf{r}-{\mathbf{r}}^{\prime }){\psi }_{n^{\prime }}\left(\mathbf{r}\right){\psi }_{m^{\prime }}\left({\mathbf{r}}^{\prime }\right).$$

(42)

The Coulomb potential itself is given by

$$V_C(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{e^2}{4\pi ϵ(|\mathbf{r}-{\mathbf{r}}^{\prime }|+\eta )},$$

(43)

where $\eta$ is a small positive regularization parameter.

Like in the continuous model, the tunneling coefficients $T_{qn}^l$ are associated to a pair of states $\{{\psi }_q^l,{\psi }_n\}$ from the lead l and the sample $S_1$. However the lattice version is much simpler:

$$T_{qn}^l=V_l{\psi }_q^{l*}\left(0_l\right){\psi }_n\left(i_l\right),$$

(44)

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

$${\psi }_q^l\left(j\right)=\frac{sin\left(q\right(j+1\left)\right)}{\sqrt{2\tau sinq}},{\epsilon }_q=2\tau cosq.$$

(45)

In the above equation $\tau$ 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 ${\epsilon }^l\left(q\right)$ scans the continuous spectrum of the semi-infinite leads ${\sigma }_l\in [-2\tau +\Delta ,2\tau +\Delta ]$ where $\Delta$ is a shift which is chosen such that ${\sigma }_l$ covers the lowest-energy many-body spectrum of the central system. The construction of the coupling coefficients $T_{qn}^l$ 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 ${\chi }_l\left(t\right)$ of the leads. It turns out that the stationary current does not depend on the last two parameters, in agreement 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 $H_{S_1}$ of the interacting system commutes with the total number operator $Q={\sum }_n{c}_n^{†}{c}_n$, its eigenstates $|\nu \rangle$ can still be labelled by the occupation $N_{\nu }$ of the non-interacting MBSs from which the state is built. Then the single index $\nu$ can be replaced by two indices, the particle number $N_{\nu }$ and an index $i_{\nu }=0,1,2,\dots$ for the ground ($i_{\nu }=0$) and excited states ($i_{\nu }>0$, where $i_{\nu }$ also counts the spin degeneracy). The notation for the interacting many-body energies is changed accordingly ${\mathcal{E}}_{\nu }\to {\mathcal{E}}_{N_{\nu }}^{\left(i_{\nu }\right)}$.

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\left(t\right)=eN\sum _{\nu ,n_{\nu }=N}\langle \nu |\rho \left(t\right)|\nu \rangle ,$$

(46)

where the sum counts all states whose total occupation $n_{\nu }=N$. Similarly one can identify the transient currents $J_{l,N}$ carried by N-particle states. These currents can be traced back form the right hand side of the GME:

$$\langle \dot{Q}\rangle =\sum _N(J_{L,N}\left(t\right)-J_{R,N}\left(t\right))=\sum _N{\dot{q}}_N\left(t\right).$$

(47)

Throughout this work we shall adopt the following sign convention for the currents associated to each lead: $J_L>0$ if the electrons flow from the left lead towards the sample and $J_R>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,\dots$). We introduce two classes of chemical potentials of the sample:

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

(48)

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

(49)

where ${\mu }_{g,N}^{\left(i\right)}$ characterizes transitions from the ground state $(N-1)$-particle configuration to various N-particle configurations. In particular ${\mu }_{g,N}^{\left(0\right)}$ describes addition processes involving ground-states with $N-1$ and N electrons while ${\mu }_{g,N}^{(i>0)}$ refers to transitions from $(N-1)$-particle ground state to excited N-particle configurations. The chemical potentials ${\mu }_{x,N}^{\left(i\right)}$ describe transitions from the 1st $(N-1)$-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 $|{\sigma }_2\rangle$. 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 $L_x=75$ nm and width $L_y=10$ nm. The lowest two spin-degenerate single-particle levels are ${\epsilon }_1=0.375$ meV and ${\epsilon }_2=3.37$ meV. The non-interacting MBSs are described by the spins of the occupied single-particle levels, e.g., $|{\sigma }_1{\overline{\sigma }}_2\rangle$ is a two-particle configuration with a spin $\sigma$ associated to the lowest single-particle state and a second electron with opposite spin orientation on the energy level ${\epsilon }_2$. Besides the usual singlet (S) and triplet (T) states we find that the Coulomb interaction induces the configuration mixing of the antiparallel configurations ${|\uparrow }_1{\downarrow }_1\rangle$ and ${|\uparrow }_2{\downarrow }_2\rangle$. More precisely, we get an interacting ground two-particle state mostly made of ${|\uparrow }_1{\downarrow }_1\rangle$ (whose weight is 0.86) and with a small component (weight 0.14) of state ${|\uparrow }_2{\downarrow }_2\rangle$. Conversely, ${|\uparrow }_1{\downarrow }_1\rangle$ is also found in the highest energy two-particle state. We stress here that the configuration mixing decreases and eventually vanishes if the gap $E_{{\uparrow }_2{\downarrow }_2}-E_{{\uparrow }_1{\downarrow }_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 ${\mu }_{g,N}^{\left(i\right)}$ lies within the bias window the corresponding state will contribute both to the transient and steady-state currents. We shall see that if ${\mu }_{g,N}^{\left(i\right)}<{\mu }_R$ the state $|i,N\rangle$ contributes only to the transient currents. Finally, when ${\mu }_{g,N}^{\left(i\right)}\gg {\mu }_L$ the state $|i,N\rangle$ is poorly populated and will not contribute to transport. Let us stress here a rather unusual transition from $|{\sigma }_2\rangle$ to the ground two-particle state which is mostly made of $|{\sigma }_1{\overline{\sigma }}_1\rangle$. The corresponding addition energy ${\mu }_{x,2}^{\left(0\right)}=2$ meV is smaller than the energy required for the usual transition $|{\sigma }_1\rangle \to |{\sigma }_1{\overline{\sigma }}_1\rangle$. This happens because of the Coulomb mixing between ${|\uparrow }_1{\downarrow }_1\rangle$ and ${|\uparrow }_2{\downarrow }_2\rangle$ which makes possible the transition from the excited single-particle state to the mixed interacting two-particle groundstate. The chemical potential ${\mu }_{x,2}^{\left(2\right)}$ describes the transition from the excited single-particle state $|{\sigma }_2\rangle$ to the triplet states.

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 ${\mu }_{g,1}^{\left(0\right)}<{\mu }_R<{\mu }_{g,1}^{\left(1\right)}<{\mu }_L<{\mu }_{g,2}^{\left(0\right)}$ (as an example we set ${\mu }_R=1$ meV and ${\mu }_L=4$ meV). Then both single-particle levels are available for tunneling but one expects that the double occupancy is excluded because ${\mu }_L<{\mu }_{g,2}^{\left(0\right)}\sim$ 5 meV. According to this scenario, more charge will accumulate on ${\epsilon }_1$, the excited states $|{\sigma }_2\rangle$ 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 ${\mu }_L<{\mu }_{g,1}^{\left(1\right)}$ as well, which suggest 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 ${\mu }_L=5.5$ meV and ${\mu }_L=4$ meV.

Figure 3 presents the evolution of the relevant populations at two values of the bias window. In Figure 3a the population $P_{1g}=P_{{\uparrow }_1}+P_{{\downarrow }_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 $|{\sigma }_2\rangle$. In particular, the ground two-particle state is populated only due to the Coulomb-induced configuration mixing.

A completely different behavior is noticed in Figure 3b. As the bias window is pushed upwards such that ${\mu }_{g,1}^{\left(1\right)}<{\mu }_R<{\mu }_{g,2}^{\left(1\right)}<{\mu }_L$ the transitions from the lowest states $|{\sigma }_1\rangle$ to two-particle states are also activated. Consequently, the population $P_{2x}$ of the excited two-particle configurations exceeds $P_{1g}$ and dominate in the steady-state regime. Note that $P_{2x}>P_{2g}$ 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_{1x}$ 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 $J_{L,N}$ and $J_{R,N}$ associated to N-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 $J_{L,1}$ and $J_{R,1}$ achieve a negative value when ${\mu }_{g,1}^{\left(0\right)}<{\mu }_{g,1}^{\left(1\right)}<{\mu }_R<{\mu }_{g,2}^{\left(0\right)}<{\mu }_L$, whereas the two-particle currents evolve to a larger positive value such that the total current $J_L$ will be positive as already shown in Figure 2b. When the bias window is shifted down to ${\mu }_L=4$ meV and ${\mu }_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 $J_L$ and that the two-particle currents display a smooth behavior. These features can be explained by looking at the charge occupations $q_N$ shown in Figure 4c,d. As ${\mu }_{g,1}^{\left(0\right)}$ and ${\mu }_{g,1}^{\left(1\right)}$ are both below ${\mu }_R=4$ meV while ${\mu }_{g,2}^{\left(0\right)}$ 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 ${\mu }_{x,2}^{\left(0\right)}$ and ${\mu }_{x,2}^{\left(1\right)}$ 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 $J_{R,2}$ slightly exceeds $J_{L,2}$.

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,64,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 $t_D=1$ meV and the time unit is expressed in units of $\hslash /t_D$. Then the currents are calculated in units of $e{t}_D/\hslash$. The tunneling rates to the four leads are all equal $V_{La}=V_{Ra}=V_{Lb}=V_{Rb}$.

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 $Q{D}_a$ is open at the initial instant $t_a=0$ and then reaches a stationary state ($J_{La}=J_{Ra}$) at some later time $T_a$. The coupling of the nearby system to its leads is switched on at $t_b>T_a$ such that the changes in the current $J_a$ 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.

As before, the interacting many-body configurations can be labeled by to the occupations of each dot according to the correspondence ${\mathcal{E}}_{\nu }\to {\mathcal{E}}_{N_{a,\nu },N_{b,\nu }}^{\left(i_{\nu }\right)}$. Here $N_{s,\nu }$ is the number of electrons in the system s associated to a many-body configuration $\nu$. If the two systems are identical the lowest chemical potentials are introduced as:

$${\mu }_g^{\left(0\right)}(N_a,N_b)={\mathcal{E}}_{N_a,N_b}^{\left(0\right)}-{\mathcal{E}}_{N_a-1,N_b}^{\left(0\right)}={\mathcal{E}}_{N_a,N_b}^{\left(0\right)}-{\mathcal{E}}_{N_a,N_b-1}^{\left(0\right)},$$

(50)

because of the degeneracy w.r.t. to the total occupation number ${\mathcal{E}}_{N_a,N_b}^{\left(0\right)}={\mathcal{E}}_{N_b,N_a}^{\left(0\right)}$ For the parameters chosen here one finds: ${\mu }_g^{\left(0\right)}(1,1)=3$ meV, ${\mu }_g^{\left(0\right)}(2,0)=4$ meV and ${\mu }_g^{\left(0\right)}(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 $J_{La}$ and $J_{Ra}$ 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 ${\mu }_{Ls}-{\mu }_{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 ${\mu }_g^{\left(0\right)}(1,1)<{\mu }_{Rs}<{\mu }_g^{\left(0\right)}(2,0)<{\mu }_{Ls}<{\mu }_g^{\left(0\right)}(2,1)$. The scenario is easy to grasp: As $Q{D}_a$ 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 ${\mathcal{E}}_{2,0}^{\left(0\right)}\leftrightarrow {\mathcal{E}}_{1,0}^{\left(0\right)}$. This behavior is observed in Figure 6a up to $t_b=150$ ps when $Q{D}_b$ is also coupled to its leads. Note also that the charge occupation of $Q{D}_a$ almost saturates at $Q_a=1.6$. As expected, for $t>t_b$ a transient current develops in $Q{D}_b$, but a simultaneous drop the $J_{La}$ and $J_{Ra}$ 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 $Q_{a,b}$ reach the same value and suggest 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 ${\mathcal{E}}_{1,1}^{\left(0\right)}$ which is not favorable for transport through any of the dots as long as ${\mu }_g^{\left(0\right)}(2,1)={\mathcal{E}}_{2,1}^{\left(0\right)}-{\mathcal{E}}_{1,1}^{\left(0\right)}$ is 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 ${\mu }_g^{\left(0\right)}(2,0)<{\mu }_{Rs}<{\mu }_g^{\left(0\right)}(2,1)<{\mu }_{Ls}$. Following the same reasoning as before one infers that now $Q{D}_a$ will enter the Coulomb blockade regime before $t=t_b$ 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 ${\mu }_g^{\left(0\right)}(2,1)={\mu }_g^{\left(0\right)}(1,2)$. This is an example of positive correlations between the two systems. Further discussions can be found in a previous publication [17].

4. Thermoelectric Transport

Until now we showed results for the charge transport driven by an electric bias of the leads due to different chemical potentials. The GME formalism allows also, in a straightforward way, the presence of a temperature bias. Instead of different chemical potentials in the left and right leads, ${\mu }_{L,R}$, one can easily consider different temperatures, $T_{L,R}$, and calculate the resulting currents after switching on the contacts between the leads on the central system. Notice that, like in the case of an electric bias, there is no requirement that the temperature bias is small, such that the nonlinear thermoelectric regime is directly accessible [69]. In addition, since the Coulomb interaction between electrons in the central system is already incorporated via the Fock space, the GME allows the inclusion of Coulomb blocking and other electron correlation effects in the thermoelectric transport [70,71].

The thermoelectric transport at nanoscale is a reach and active topic within the context of the modern quantum thermodynamics, partly motivated by novel ideas on the conversion of wasted heat into electricity, and partly by the characterization of nanoscale system by methods complementary to pure electric transport [72]. For example, an effect specific to nanosystems is the sign change of the thermoelectric current or voltage when the electronic energy spectrum consists of discrete levels. This effect was predicted in the early 90’ [73] and detected experimentally for quantum dots [74,75,76] and molecules [77]. This means that thermoelectric current in a nanoelectronic system may flow from the hotter contact to the colder one, but also from the colder to the hotter, although the second possibility might appear counter-intuitive. A simple explanation of this sign change of the current is that in a nanoscale system with discrete resonances the current can be seen as having two components, one carried by populated states above the Fermi energy, and another one carried by depopulated states below it. By analogy with a semiconductor, the former states correspond to electrons in the conduction band and the later states to holes in the valence band. Whereas an electric bias drives the electric currents due to particles and holes in the same direction, such that they always add up, a thermal bias drives them in opposite directions, such that the net current is their difference, which can be positive, negative, or zero.

We can describe this effect with the GME, first assuming a simple model with unidimensional and discretized leads, and just a single site in between them as central system. By using the Markov approximation one can show analytically that the current in the leads, in the steady state, are obtained as [71]

$$J_{L,R}=\frac{1}{{\tau }^2}\frac{V_L^2{V}_R^2}{V_L^2+V_R^2}\left[f_L\left(E\right)-f_R\left(E\right)\right],$$

(51)

where $V_{L,R}$ are the coupling parameters of the leads with the central site, $\tau$ is the hopping energy on the leads, and E is the energy of the central site. We see that the sign of the current depends on the difference between the Fermi energies in the leads at the resonance energy,

$$f_l\left(E\right)=\frac{1}{e^{(E-{\mu }_l)/k_B{T}_l}+1},l=L,R.$$

(52)

Thus, in the presence of a thermal bias, say $T_L>T_R$, but in the absence of an electric bias, i.e., ${\mu }_L={\mu }_R$, the current is zero and changes sign around ${\mu }_l=E$. In addition, the current may also vanish if the chemical potential in the leads is sufficiently far from the resonance such that the two Fermi functions are both close to zero or one. Which means that if the central system has more resonant energies the current may also change sign when ${\mu }_l$ is somewhere between two of them.

In Figure 7 we show an example of thermoelectric currents calculated with the GME, using the same model as in Section 3. The lowest single-particle levels having energies ${\epsilon }_1=0.375$ meV and ${\epsilon }_2=3.37$ meV are followed by the two-particle singlet state with $E_s=5.39$ meV and triplet with $E_t=5.62$ meV, and then by another excited two-body state with zero spin with energy $E_x=10.5$ meV. We consider temperatures $k_B{T}_L=0.5$ meV and $k_B{T}_R=0.05$ meV in the left and right lead, respectively (or $T_L=5.8$ K and $T_R=0.58$ K), and equal chemical potentials. In Figure 7a one can see the time dependence of the currents in the leads after they are coupled to the central system, for two values of the chemical potentials, 4.8 meV and 5.4 meV, selected on each side of the singlet state. Compared to the results shown in Section 3 here we increased the coupling parameters between the leads and the central system 1.4 times, such that the steady state is reached sooner.

As predicted by Equation (52), the currents in the steady state have opposite sign. But in fact, as shown by the red curve of Figure 7b, here we do not resolve the energy interval between the singlet and triplet states with $k_B{T}_L>E_t-E_s=0.23$ meV, such that we obtain one single (common) sign change for these two levels (or “resonances”). Next, by increasing the chemical potential within the larger gap between $E_s$ and $E_x$ the current in the steady state approaches zero and changes sign again, for ${\mu }_l\approx 7.0$ meV, and for ${\mu }_l\approx 7.8$ meV when the temperature of the hot lead is doubled, $T_L=11.5$ K.

By varying the chemical potential below the singlet energy $E_s$ we obtain a similar decreasing trend of the current, except that now there is no sign change close to the energy ${\epsilon }_2=3.37$ meV, but only a succession of minima and maxima. The reason is the level broadening due to the coupling of the central system with the leads [71]. Still, from such data one can observe experimentally the charging energy, as the interval between consecutive maxima, or minima, or mid points between them [77].

In the present review we show only the thermoelectric current, which corresponds to a short-circuit experimental setup, i.e., a circuit without a load. To obtain a voltage with the GME method one has to simulate a load by considering also a chemical potential bias. Thus, one can obtain the open-circuit voltage, which corresponds to that electric bias ${\mu }_R-{\mu }_L$ which totally suppresses the thermoelectric current, or the complete I-V characteristic of the “thermoelectric device”. Interestingly, the sign change of the thermoelectric current or voltage can also be obtained by increasing the temperature of the hot lead, while keeping the other lead as cold as possible [70,78,79,80].

A novel example of sign reversal of the thermoelectric current has been recently predicted in tubular nanowires, either with a core-shell structure or made of a topological insulator material, in the presence of a transversal magnetic field [81]. In this case the energy spectra are continuous, but organized in subbands which are nonmonotonic functions of the wavevector along the nanowire, yielding a transmission function nonmonotonic with the energy, and the reversal of the thermoelectric current, even in the presence of moderate perturbations [82,83].

5. Electron Transport through Photon Cavities

5.1. The Electron-Photon Coupling

From the beginning our effort to model electron transport through a nano scale system placed in a photon cavity has been geared towards systems based on a two-dimensional electron gas in GaAs or similar heterostructures. We have emphasized intersubband transitions in the conduction band, active in the teraherz range, in anticipation of experiments in this promising system [84].

Here, subsystem $S_1$ is a two-dimensional electronic nanostructure placed in a static (classical) external magnetic field. The leads are subjected to the same homogeneous external field. The electronic nanostructure, via split-gate configuration, is parabolically confined in the y-direction with a characteristic frequency ${\Omega }_0$. The ends of the nanostructure in the x-direction at $x=\pm L_x/2$ are etched, forming a hard-wall confinement of length $L_x$. The external classical magnetic field is given by $\mathbf{B}=B\widehat{\mathbf{z}}$ with a vector potential $\mathbf{A}=(-By,0,0)$. The single-particle Hamiltonian reads:

$$\begin{array}{ccc}\hfill {\widehat{h}}_{S_1}^{\left(0\right)}& =& \frac{1}{2m}{\left(\mathbf{p}+q\mathbf{A}\right)}^2+\frac{1}{2}m{\Omega }_0^2{y}^2\hfill \\ \hfill & =& \frac{1}{2m}{p}_x^2+\frac{1}{2m}{p}_y^2+\frac{1}{2}m{\Omega }_w^2{y}^2+i{\omega }_{c}y{p}_x,\hfill \end{array}$$

(53)

where m is the effective mass of an electron, $-q$ its charge, $\mathbf{p}$ the canonical momentum operator, ${\omega }_c=qB/m$ is the cyclotron frequency and ${\Omega }_w=\sqrt{{\omega }_c^2+{\Omega }_0^2}$ is the modified parabolic confinement. The spin degree of freedom is included with either a Zeeman term added to the Hamiltonian [85], or with Rashba and Dresselhaus spin orbit interactions, additionally [86].

$H_{S_2}$ is simply the free field photon term for one cavity mode and by ignoring the zero point energy can be written as $H_{S_2}=\hslash {\omega }_p{a}^{†}a$ where $\hslash {\omega }_p$ is the single photon energy and a ($a^{†}$) is the bosonic annihilation (creation) operator. The electron-photon interaction term $V_{\mathrm{el}-\mathrm{ph}}$ can be split into two terms $V_{\mathrm{el}-\mathrm{ph}}=V_{\mathrm{el}-\mathrm{ph}}^{\left(1\right)}+V_{\mathrm{el}-\mathrm{ph}}^{\left(2\right)}$ where

$$\begin{array}{ccc}\hfill V_{\mathrm{el}-\mathrm{ph}}^{\left(1\right)}& :=& {\displaystyle \sum _{ij}}{\displaystyle \sum _{i,j}}〈{\psi }_i\left|\frac{q}{2m}\left(\pi ·{\mathbf{A}}_{\mathrm{EM}}+{\mathbf{A}}_{\mathrm{EM}}·\pi \right)\right|{\psi }_j〉c_i^{†}{c}_j\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill V_{\mathrm{el}-\mathrm{ph}}^{\left(2\right)}& :=& {\displaystyle \sum _{ij}}{\displaystyle \sum _{i,j}}〈{\psi }_i\left|\frac{q^2}{2m}{\mathbf{A}}_{\mathrm{EM}}^2\right|{\psi }_j〉c_i^{†}{c}_j,\hfill \end{array}$$

(55)

with $\pi \equiv \mathbf{p}+q\mathbf{A}$ the mechanical momentum. The term in Equation (54) is the paramagnetic interaction, whereas the diamagnetic term is defined by Equation (55). By assuming that the photon wavelength is much larger than characteristic length scales of the system one can approximate the vector potential amplitude to be constant over the electronic system. Let us stress here that, in contrast to the usual dipole approximation, we will not omit the diamagnetic electron-photon interaction term. Then the vector potential is written as:

$${\mathbf{A}}_{\mathrm{EM}}\simeq \widehat{\mathbf{e}}{A}_{\mathrm{EM}}\left(a+a^{†}\right)=\widehat{\mathbf{e}}\frac{{\mathcal{E}}_c}{q{\Omega }_w{a}_w}\left(a+a^{†}\right),$$

(56)

where $\widehat{\mathbf{e}}$ is the unit polarization vector and ${\mathcal{E}}_c\equiv q{A}_{\mathrm{EM}}{\Omega }_w{a}_w$ is the electron-photon coupling strength. For a 3D rectangular Fabry Perot cavity we have $A_{\mathrm{EM}}=\sqrt{\hslash /\left(2{\omega }_{p}V{ϵ}_0\right)}$ where V is the cavity volume. Linear polarization in the x-direction is achieved for a ${\mathbf{TE}}_{011}$ mode, and in the y-direction with a ${\mathbf{TE}}_{101}$ mode.

Using the approximation in Equation (56), the expressions for the electron-photon interaction in Equations (54) and (55) are greatly simplified by pulling ${\mathbf{A}}_{\mathrm{EM}}$ in front of the integrals. For the paramagnetic term, we get

$$V_{\mathrm{el}-\mathrm{ph}}^{\left(1\right)}\simeq {\mathcal{E}}_c\left(a+a^{†}\right)\sum _{ij}{g}_{ij}{c}_i^{†}{c}_j.$$

(57)

where we introduced the dimensionless coupling between the electrons and the cavity mode

$$g_{ij}=\frac{a_w}{2\hslash }\widehat{\mathbf{e}}·\int d\mathbf{r}\left[{\psi }_i^{*}\left(\mathbf{r}\right)\left\{\pi {\psi }_j\left(\mathbf{r}\right)\right\}+\left\{\pi {\psi }_i^{*}\left(\mathbf{r}\right)\right\}{\psi }_j\left(\mathbf{r}\right)\right].$$

(58)

As for the diamagnetic term, we get

$$V_{\mathrm{el}-\mathrm{ph}}^{\left(2\right)}\simeq \frac{{\mathcal{E}}_c^2}{\hslash {\Omega }_w}\left[\left(a^{†}a+\frac{1}{2}\right)+\frac{1}{2}\left(a^{†}{a}^{†}+aa\right)\right]{\mathcal{N}}^e,$$

(59)

where ${\mathcal{N}}^e$ is the number operator in the electron Fock space. Note that $V_{\mathrm{el}-\mathrm{ph}}^{\left(2\right)}$ does not depend on the photon polarization or geometry of the system in this approximation. We do not use the rotating wave approximation as in our multilevel systems even though a particular electron transition could be in resonance with the photon field we want to include the contribution form others not in resonance.

For the numerical diagonalization of $H_S$ we shall use the lowest $N_{\mathrm{mesT}}\ll N_{\mathrm{mes}}$ IMBS of $H_{S_1}$ and photon states containing up to $N_{\mathrm{EM}}$ photons, resulting in a total of $N_{\mathrm{mesT}}\times (N_{\mathrm{EM}}+1)$ states in the ‘free’ basis $\{|\nu ,j\rangle \}$.

5.2. Results

Groups modeling the near resonance interaction of one cavity mode with a two level electronic system have expressed the importance of using a large enough, or the correct type, of a photon basis in the strongly interacting regime [87,88]. In many level systems where wavefunction and geometric effects are accounted for our experience is that convergence in numerical diagonalization is more sensitive to proper truncation of the electronic sector of the Fock many-body space. This reflects the polarizability of the electric charge by a cavity field in the construction of the photon-dressed electronic states. At the same time the inclusion of the diamagnetic interaction curbs the need for states with a very high photon number [42,50,89].

The polarizability of the first photon replica of the two-electron ground state is displayed in Figure 8 as a function of $g_{\mathbf{EM}}$, the photon energy $\hslash \omega$ and its polarization [50].

The polarizability is nonlinear, anisotropic, and largest for the cavity photon close to a resonance with the confinement energy in the y-direction.

A Rabi oscillation of two electrons in the double quantum dot system embedded in the short quantum wire leads to oscillating charge with time in the system. The oscillating probability of charge presence in the contact areas of the short wire thus lead to oscillations in the current leaving the system through the left and right leads [33], see Figure 9.

Alternatively, one may view this as the consequence of the Rabi resonance entangling two states with different tunneling probability to the leads.

In the transient or the late transient regime we have used the non-Markovian GME to investigate several results: Thorsten Arnold et al. used a time-convolution-less (TCL) version of the GME to study the effects of magnetic field and photons [46] on the transport of interacting electrons through a quantum ring with spin-orbit interactions in a photon cavity with circular [86] and linear polarization [90]. Aharonov-Bohm oscillations were established in the time-dependent transport through a ring structure with additional vortexes in the contact region of the quantum wire. x-polarized photons with energy 0.3 meV attenuate the Aharonov-Bohm oscillations over a broad range of magnetic field, but y-polarized photons influence the transport in a more complex fashion. The oscillations are generally attenuated, but one oscillation peak is split and the charge current is enhanced at a magnetic field corresponding to a half-integer flux quantum [46]. With the spin-orbit interactions the spin polarization and the spin photo currents of the quantum ring are largest for circularly polarized photon field and a destructive Aharonov–Casher (AC) phase interference. The dip in the charge current caused by the destructive AC phase becomes threefold under the circularly polarized photon field as the interaction of the angular momentum of the electron and the spin angular momentum of the light create a many-body level splitting [86]. The detailed balance between the para- and the diamagnetic electron-photon interactions has been studied for an electron in the quantum ring structure when excited by a short classical dipole pulse [47].

Nzar Rauf Abdullah et al. have used the GME formalism to investigate photon assisted transport [91], photon mediated switching in nanostructures [48,49,92], the balancing of magnetic and forces caused by cavity photons [93], cavity-photon affected thermal transport [94,95], and the influence of cavity photons on thermal spin currents in a system with spin orbit interactions [96,97].

6. Steady-State

The investigation of the time dependent transport of electrons through a photon cavity soon made it clear that for the continuous model the inherent time scales can lead to relaxation times far beyond what is accessible with simple integration of the GME [32,45,46,49,91]. The underlying cause for the diverse relaxation times is on one hand electron tunneling rates affected by the shape or geometry of the system and the condition of weak coupling. Different many-body states can have a high or low probability for electrons to be found in the contact areas of the central system. On the other hand are slow rates of FIR or terahetz active transitions, that are furthermore affected by the geometry of the wavefunctions of the corresponding final and initial states. In addition, the cavity decay, or coupling to the environment, affects relaxation times as we address below [98]. To avoid confusion it is important to remember that we calculate the eigenstates of the closed central system, the interacting electron and photon system, and the opening up of the system to the leads or the external photon reseroir is always a neccessary triggering mechanism for all transitions later in time, photon active or not.

6.1. The Steady-State Limit

In order to investigate the long-time evolution and the steady state of the central system under the influences of the reservoirs we resort to a Markovian version of the GME, whereby we assume memory effects in the kernel of the GME (26) to vanish, relinquishing the reduced density operator local in time enabling the approximation [99]

$${\int }_0^{\infty }ds{e}^{is(E_{\beta }-E_{\alpha }-ϵ)}\approx \pi \delta (E_{\beta }-E_{\alpha }-ϵ),$$

(60)

where a small imaginary principle part is ignored. We have furthermore assumed instant lead-system coupling at $t=0$ with ${\chi }_l\left(t\right)=\theta \left(t\right)$, the Heaviside unit step function, in Equation (8) for $H_T$. In order to transform the resulting Markovian equation into a simpler form we use the vectorization operation [100], that stacks the columns of a matrix into a vector, and its property

$$\mathrm{vec}\left(A\rho B\right)=(B^T\otimes A)\mathrm{vec}\left(\rho \right)$$

(61)

through which the reduced density matrix can always be moved to the right side of the corresponding term, and a Kronecker product has been introduced with the property $B\otimes A=\left\{B_{\alpha ,\beta }A\right\}$. The Kronecker product of two $N_{\mathrm{mes}}\times N_{\mathrm{mes}}$ matrices results in a $N_{\mathrm{mes}}^2\times N_{\mathrm{mes}}^2$ matrix, and effectively the vectorization has brought forth that the natural space for the Liouville-von Neumann equation is not the standard Fock space of many-body states, but the larger Liouville space of transitions [101,102,103].

No further approximations are used to attain the Markovian master equation and due to the complex structure of the non-Markovian GME we have devised a general recipe published elsewhere [99] to facilitate the analytical construction and the numerical implementation. The Markovian master equation has the form

$${\partial }_t\mathrm{vec}\left(\rho \right)=\mathcal{L}\mathrm{vec}\left(\rho \right),$$

(62)

and as the non-Hermitian Liouville operator $\mathcal{L}$ is independent of time the analytical solution of Equation (62) can be written as

$$\mathrm{vec}\left(\rho \left(t\right)\right)=\left\{\mathcal{U}[exp\left({\mathcal{L}}_{\mathrm{diag}}t\right)]\mathcal{V}\right\}\mathrm{vec}\left(\rho \left(0\right)\right),$$

(63)

in terms of the matrices of the column stacked left $\mathcal{U}$, and the right $\mathcal{V}$ eigenvectors of $\mathcal{L}$

$$\mathcal{L}\mathcal{V}=\mathcal{V}{\mathcal{L}}_{\mathrm{diag}},\mathrm{and}\mathcal{U}\mathcal{L}={\mathcal{L}}_{\mathrm{diag}}\mathcal{U},$$

(64)

obeying

$$\mathcal{U}\mathcal{V}=I,\mathrm{and}\mathcal{V}\mathcal{U}=I.$$

(65)

Special care has to be taken in the numerical implementation of this solution procedure as many software packages use another normalization for $\mathcal{U}$ and $\mathcal{V}$. Calculations in the Liouville space using (63) are memory (RAM) intensive, but bring several benefits: No time integration combined with iterations is needed, thus time points can be selected with other criteria in mind. The solution is thus convenient for long-time evolution, that is not easily accessible with numerical integration. The complicated structure of the left and right eigenvector matrices for a complex system with nontrivial geometry makes Equation (63) the best choice to find the properties of the steady state by monitoring the limit of it as time gets very large. The complex eigenvalue spectrum of the Liouville operator $\mathcal{L}$ reveals information about the mean lifetime and energy of all active transitions in the open system, and the zero eigenvalue defines the steady state.

In the steady state the properties of the system do not change with time, but underneath the “quiet surface” many transitions can be active to maintain it. The best experimental probes to gauge the underlying processes are measurements of noise spectra for a particular physical variable. They are available through the two-time correlation functions of the respective measurable quantities. For a Markovian central system weakly coupled to reservoirs the two-time correlation functions can be calculated applying the Quantum Regression Theorem (QRT) [104,105] stating that the the equation of motion for a two-time correlation function has the same form as the Markovian master equation (62) for the operator [106]

$$\chi \left(\tau \right)={\mathrm{Tr}}_{\mathrm{R}}\left\{e^{-iH\tau /\hslash }X{\rho }_{\mathrm{T}}\left(0\right)e^{+iH\tau /\hslash }\right\},$$

(66)

where H is the total Hamiltonian of the system, ${\rho }_{\mathrm{T}}$ its density operator, and the trace is taken with respect of all reservoirs. For photon correlations $X=a+a^{†}$ as in [50], or $X=Q{\Lambda }^l$ for current correlations as in [107], where $Q={\sum }_i{c}_i^{†}{c}_i$ is the fermionic charge operator and ${\Lambda }^l$ is the Liouville dissipation operator for lead l. The structure of $\chi$ (66) indicates that the two-time correlation function is then

$$\langle X\left(\tau \right)X\left(0\right)\rangle ={\mathrm{Tr}}_{\mathrm{S}}\left\{X\left(0\right)\chi \left(\tau \right)\right\},$$

(67)

with

$$\mathrm{vec}\left(\chi \left(\tau \right)\right)=\left\{\mathcal{U}[exp\left({\mathcal{L}}_{\mathrm{diag}}t\right)]\mathcal{V}\right\}\mathrm{vec}\left(\chi \left(0\right)\right).$$

(68)

The left side of Equation (67), the two-time correlation function, is written in the Heisenberg picture, in contrast to the Schrödinger picture used elsewhere in the article. The Fourier spectral density for the photon two-time correlation function is denoted by $S\left(E\right)$, and for the current-current correlation the corresponding Fourier spectral density denoted by $D_{l{l}^{\prime }}\left(E\right)$, where l and $l^{\prime }$ refer to L and R, the Left and Right leads.

6.2. Results

To date we have used the Markovian version of the master equation to investigate properties of the steady state, and how the system with electrons being transported through a photon cavity reaches it. We assume GaAs parameters with effective mass $m=0.067m_e$, effective relative dielectric constant ${ϵ}_r=12.3$, and effective Landé g-factor $g=-0.44$. The characteristic energy of the parabolic confinement of the semi-infinite leads and the central system in the y-direction is $\hslash {\Omega }_0=2.0$ meV. The length of the short quantum wire is $L_x$, and the overall coupling coefficient for the leads to the system is $g_{LR}{a}_w^{3/2}=0.124$ meV.

We start with a central system made of a finite parabolic quantum wire without any embedded quantum dots. Figure 10 demonstrates that the approach to build and solve the Markovian master Equations (62)–(63) works for an interacting system with 120 many-body states participating in the transport [108].

The upper right panel displays the properties of the lowest 32 many-body states at the plunger gate voltage $V_g=-1.6$ mV. With ${\mu }_L=1.4$ meV and ${\mu }_L=1.1$ meV there are 8 states below the bias window and five states within it. In the bias window is one spin singlet two-electron state (the two-electron ground state) and two spin components of two one-electron states with a non-integer mean photon content indicating a Rabi splitting. The upper left panel of Figure 10 show the mean electron and photon numbers in the central system when it is initially empty. With a very low coupling, $g_{\mathrm{EM}}=1\times {10}^{-6}$ meV, between the electrons and photons, the charging is very slow with the probability approaching unity around $t\approx {10}^8$ s. With increasing $g_{\mathrm{EM}}$ the charging becomes faster, and during the phase the mean photon number in the system rises. The lower panels of Figure 10 reveal what is happening. With the low photon coupling (lower left panel) electrons tunnel non-resonantly into the two spin components of the ground state, $|1)$ and $|2)$ as the vacuum state $|3)$ looses occupation, and to a small fraction the two-electron state $|9)$ gets occupied. When the coupling of the electrons and the photons is not vanishingly small (lower right panel) the charging of the system takes a different rout. The finite $g_{\mathrm{EM}}$ allows the incoming electron to enter the Rabi-split one-electron states in the bias window as these are a linear combination of electron states with a different photon number. This explains the growing mean number of photons in the system for intermediate times. These states are eigenstates of the central system, but not of the open system, so at a later time they decay into the the one- and two-electron ground states as before bringing the system into the same Coulomb blocked steady state as before. We thus observe electromagnetically active transitions in the system in an intermediate time regime [108].

The on-set of the steady state regime is difficult to judge only from the shape of the charge being accumulated in the system or the current through or into it as a function of time [85]. For a system of two parallel quantum dots embedded in a short quantum wire ($L_x=150$) nm the charging and the current as functions of time look the same (see Figures 4 and 5 in ref. [85]), but when the occupation of the eigenstates of the closed system is analyzed, see Figure 11, a clear difference is seen for the approach to the steady state depending on whether the initial state contains only one or no photon [85]. In the case of neither photon nor an electron in the cavity initially an electron tunnels into the system into the two spin components of the one-electron ground state, which happens to be in the bias window for $V_g=-2.0$ mV. Thus, the steady state is a combination of the empty state and these two one-electron states. In the case of one photon and no electron initially in the system an electron tunnels non-resonantly into the 1-electron states $|8)$ and $|9)$ with energy slightly below 2 meV, and thus well above the bias window. The mean photon content of these states is close to unity and at a later time the electron ends up in the two spin components of the one-electron ground state via a radiative transition [85]).

Note that the “irregularly” looking structure around $t\approx 2000$ ps will be addressed below. Please note that the numbering of interacting many-body state depends on the structure of the system, and the plunger gate voltage $V_g$.

In the steady state all the mean values of the open system have reached a constant value. In order to query about the active underlying processes it is necessary to calculate the spectral densities of the photon or current correlations. We present these for the central system consisting of a short quantum wire ($L_x=150$) nm with two embedded quantum dots in Figure 12 (see refs. [107,109]).

Importantly we show in Ref. [50] that both the paramagnetic and the diamagnetic electron-photon interactions can lead to a Rabi resonance. The resonance for the diamagnetic interactions is much smaller, but the symmetry of the two parallel quantum dots leads to selection rules where for x-polarized cavity photon field the paramagnetic interaction is blocked, but both are present for the y-polarized field. Here, the active states are the one-electron ground state and the first excited one-electron state, with which the first photon replica of the ground state interacts for $\hslash \omega =0.72$ meV. The spectral density of the photon-photon two-time correlation function, $S\left(E\right)$ seen in the left panel of Figure 12 shows one peak at the energy of the cavity mode $\hslash \omega =0.72$ meV, and two side peaks for the y-polarization. The central peak is the ground state state electroluminescence and the side peaks are caused by the Rabi-split states [34,109,110]. Here, we observe the ground state electroluminescence even though the electron-photon coupling is not in the ultra strong regime, as we diagonalize the Hamiltonian in a large many-body Fock space instead of applying conventional perturbative calculations.

For the x-polarized cavity field we find a much weaker ground state electroluminescence caused by the diamagnetic electron-photon interaction [109]. In addition, we identify these effects for the fully interacting two-electron ground state, where they are partially masked by many concurrently active transitions. The spectral density for the current–current correlation functions $D_{l{l}^{\prime }}\left(E\right)$ displayed in the right panel of Figure 12 show only peaks at the Rabi-satellites, as could be expected [107]. An inspection of $D_{l{l}^{\prime }}\left(E\right)$ over a larger range of energy reveals more transitions active in maintaining the steady state, both radiative transitions and non radiative [107]. Moreover, we notice that when the steady state is not in a Coulomb blocking regime the spectral density of the current-current correlations always shows a background to the peaks with a structure reminiscent of a $1/f$ behavior, that is known in multiscale systems.

An “irregularly” looking structure in the mean occupation, the current current, and the mean number of electrons and photons. This is a general structure seen in all types of central system we have investigated in the continuous model. In Figure 13 we analyze it in a short parabolically confined quantum wire of length $L_x=180$ nm with two asymmetrically embedded quantum dots [111].

An increased number of time points on the logarithmic scale shows regular oscillations. A careful analysis reveals two independent oscillations: A spatial charge oscillation between the quantum dots with the Rabi frequency in the system, and a still slower nonequilibrium oscillation of the spin populations residing as the system is brought to a steady state [111].

The steady-state Markovian formalism has been used to investigate oscillations in the transport current as the photon energy or the electron–photon coupling strenght are varied with or without flow of photons from the external reservoir [112,113]. Moreover, the formalism has been used to establish the signs of the Purcell effect [114] in the transport current [98].

In light of the experimental interest of using a two-dimensional electron gas in a GaAs heterostructure [84] we have calculated the exact matrix elements for the electron-photon interaction taking into account the spatial variation of the vector field $\mathbf{A}$ of the electronic system. This is a small correction in most cases but may be important when studying high order transitions or nonperturbational effects caused by the photon field. This has led us to discover a very slow high order transition between the ground states of two slightly dissimilar quantum dots [115].

The fist steps have been taken to investigate thermoelectric effects in the central system coupled to cavity photons, in the steady state. In a short quantum wire with one embedded quantum dot in the resonant regime, an inversion of thermoelectric current is found caused by the Rabi-splitting. The photon field can change both the magnitude and the sign of the thermoelectric current induced by the temperature gradient in the absence of a voltage bias between the leads [116].

7. Summary

It goes without saying that as transport experiments at nanoscale become more involved the formal tools must be suitably extended or adapted. In particular, the unavoidable charging and correlation effects at finite bias pushed the theoretical calculations from the very convenient single-particle (or at most mean-field) Landauer–Büttiker picture to the complicated many-body perturbation theory of the non-equilibrium Keldysh–Green’s functions [117].

Here we summarized some results on time-dependent transport in open interacting systems which argue for the similar idea that if one looks for transient effects and dynamics of excited states the simple rate equation approach must be extended to the non-markovian generalized master equation.

The GME we used in all examples is constructed and solved w.r.t the exact many-body states of the central open system and can be therefore implemented numerically without major changes to study both Coulomb-interacting and hybrid systems where the fermion-boson interaction is crucial, like QD-cavity systems or nano-electromechanical systems. A consistent derivation of the GME the full knowledge of the eigenvalues and eigenfunctions of complicated interacting Hamiltonian (e.g., cavity-coupled systems must be described by ‘dressed’ states). With very few exceptions coming from quantum optics (i.e., the Jaynes-Cummings model for two-level or $\Lambda$ and V three-level systems) such a task can only be achieved via numerically exact diagonalization of large matrices, especially for elecron-photon systems. To bypass this difficulty we proposed and succesfully used a stepwise diagonalization procedure.

The dynamics of excited states in a quantum wire, the onset of current-current correlations for a pair of electrostatically coupled quantum dots and thermoelectric effects were presented within a simple lattice model which however captures the relevant physics.

When turning to QED-cavity system we developed the GME within a continuous model which accounts for the geometrical details of the sample and of the contact regions. Moreover, the calculations were performed by taking into account both the paramagnetic and diamagnetic contributions to the electron-photon coupling and without relying on the rotating-wave approximation. This is an important step beyond the Jaynes–Cummings model. Also, the number of many-body stated needed in the calculations increased considerably. Thus, the accuracy of the stepwise numerical diagonalization had to be carefully discussed. Finally, for systems with long relaxation time a markovian version of GME was proposed and implemented via a clever vectorization procedure.

We end this review by pointing out possible improvements of the GME method and some of its future applications. At the formal level, perhaps the most challenging upgrade is the inclusion of time-dependent potentials describing laser pulses or microwave driving signals. Provided this is succesfully achieved, one could study transport through driven nano-electromechanical systems (NEMS) or the physics of Floquet states emerging in strongly driven systems [118,119]. Let us mention here that at least for closed systems (i.e., not connected to particle reservoirs) studies based on Floquet master equations for two-level system are already available [120,121]. As for more immediate applications we aim at the theoretical modeling of transport in Tavis–Cummings systems, motivated by the recent observation of state readout in a system of distant coupled quantum dots individually connected to a pair of leads and interacting via cavity photons [29].