**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 *HSR* 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 *nr* particle reservoirs characterized by chemical potentials and temperatures {*μl*, *Tl*}, *l* = 1, 2, ..., *nr*, 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 {*<sup>ω</sup>k*}. Let F*S*1 and F*S*2 be the Fock spaces associated to the two systems. Typically F*S*1 is a set of interacting many-body configurations of the electronic system whereas F*S*2is 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\_{\mathbb{S}} = H\_{\mathbb{S}\_1} + H\_{\mathbb{S}\_2} + V. \tag{1}$$

In this work *HS*1 will describe various Coulomb-interacting structures: a single quantum dot, a 2D wire or parallel quantum dots. We shall denote by |*ν* and *Eν* the many-body configurations and eigenvalues of *HS*1 , that is one has *HS*1 |*ν* = *<sup>E</sup>ν*|*ν*. *HS*1 can be equally expressed in terms of creation and annihilation operators {*c*† *nσ*, *cnσ*} associated to a spin-dependent single-particle basis {*ψnσ*} of a single-particle Hamiltonian *h* (0) *S*1(see the next sections for specific models), such that:

$$H\_{\mathbb{S}\_1} = H\_{\mathbb{S}\_1}^{(0)} + \mathcal{W}\_{\prime} \tag{2}$$

where *H*(0) *S*1 is the 2nd quantized form of *h* (0) *S*1 and *W* is the Coulomb interaction. Similarly, the eigenstates and eigenvalues of the second subsystem *S*2 will be denoted by |*j* and *ej* such that *HS*2 |*j* = *ej*|*j*. 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(t) = H\_S + H\_R + H\_{SR}(t) \; := H\_0 + H\_{SR}(t), \tag{3}$$

where the system-reservoir coupling *HSR* collects the coupling to the leads (*HT*) and the coupling of a bosonic mode to a thermal or leaky bosonic environments (*HE*):

$$H\_{SR}(t) = H\_T(t) + H\_E. \tag{4}$$

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

$$H\_R = H\_{\text{leads}} + H\_{\text{both}} \tag{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 *HS* 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 *ε<sup>l</sup>*(*q*) and *<sup>ψ</sup>lqσ* 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 *<sup>c</sup>*†*qlσ*/*cql<sup>σ</sup>* associated to the single particle states, we can write: 

$$H\_{\rm leads} = \sum\_{l} H\_{l} = \int dq \sum\_{\sigma} \epsilon^{l}(q) c\_{qlr}^{\dagger} c\_{qlr} \,. \tag{6}$$

As for the bosonic bath, it is described by a collection of harmonic oscillators with frequencies *ωk* and by corresponding creation/annihilation operators *<sup>b</sup>*†*k*/*bk*:

$$H\_{\rm bath} = \sum\_{k} \hbar \omega\_{k} b\_{k}^{\dagger} b\_{k}. \tag{7}$$

The tunneling Hamiltonian has the usual form:

$$H\_T(t) = \sum\_{l} \sum\_{nr} \int dq \chi\_l(t) \left( T\_{qn}^l c\_{qlr}^\dagger c\_{nr} + h.c \right), \tag{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 *σ* will be henceforth tacitly merged with the single-particle index *n* and restored if needed.

The time-dependent switching functions *<sup>χ</sup>l*(*t*) 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 *<sup>χ</sup>l*(*t*) = *θ*(*t* − *t*0) where *<sup>θ</sup>*(*x*) 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 *Tlqn* describes the tunneling strength between a state with momentum *q* of the lead *l* and the state *n* of the isolated sample with wavefunctions *ψ<sup>n</sup>*. 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 W of the open system obeys the Liouville–von Neumann equation:

$$i\hbar \frac{\partial \mathcal{W}(t)}{\partial t} = \mathcal{L}(t)\mathcal{W}(t), \quad \mathcal{W}(t\_0) = \rho\_S(t\_0) \odot \rho\_{\mathcal{R}} \tag{9}$$

where:

$$
\mathcal{L}(t) = \mathcal{L}\_0 + \mathcal{L}\_{SR}, \quad \mathcal{L}\_0 \cdot = [\mathcal{H}\_{0\prime} \cdot]. \tag{10}
$$

We also introduce the notations:

$$
\mathcal{L}\_{\mathbb{S}^\cdot} = [H\_{\mathbb{S}\prime} \cdot], \quad \mathcal{L}\_{SR} \cdot = [H\_{SR\prime} \cdot]. \tag{11}
$$

The Nakajima–Zwanzig projection formalism leads to an equation of motion for the reduced density operator *ρ*(*t*) = Tr*R*{W}. The initial state W0 := W(*<sup>t</sup>*0) factorizes as:

$$\mathcal{W}\_0 = \rho\_0 \otimes \rho\_{\text{leads}} \otimes \rho\_{\text{both}} := \rho\_0 \otimes \rho\_{\mathcal{K}} \tag{12}$$

where the equilibrium density operator of the leads reads:

$$\rho\_{\text{leads}} = \prod\_{l} \frac{e^{-\beta\_l (H\_l - \mu\_l N\_l)}}{\text{Tr}\_l \left\{ e^{-\beta\_l (H\_l - \mu\_l N\_l)} \right\}},\tag{13}$$

and *βl* = 1/*kBTl*, *μl* and *Nl* denote the inverse temperature, chemical potential and the occupation number operator of the lead *l*. Similarly,

$$\rho\_{\rm bath} = \prod\_{k} e^{-\hbar \omega\_{k} b\_{k}^{\dagger} b\_{k} / k\_{B}T} (1 - e^{-\hbar \omega\_{k} / k\_{B}T}) . \tag{14}$$

Finally, *ρ*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 \cdot = \rho\_R \text{Tr}\_R\{\cdot\}, \quad Q = 1 - P. \tag{15}$$

It is straightforward to check the following properties:

$$P\mathcal{L}\_S \quad = \quad \mathcal{L}\_S P, \quad P\mathcal{L}\_{SR}(t)P = 0. \tag{16}$$

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

$$i\hbar \mathcal{P} \frac{\partial \mathcal{W}(t)}{\partial t} = -P \mathcal{L}(t) P \mathcal{W}(t) + P \mathcal{L}(t) Q \mathcal{W}(t) \tag{17}$$

$$i\hbar \dot{Q} \frac{\partial \mathcal{W}(t)}{\partial t} = -Q \mathcal{L}(t) Q \mathcal{W}(t) + Q \mathcal{L}(t) P \mathcal{W}(t),\tag{18}$$

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

$$Q\mathcal{W}(t) = \frac{1}{i\hbar} \int\_{t\_0}^{t} ds \, T \exp\left\{-\frac{i}{\hbar} \int\_{s}^{t} ds' \mathcal{QL}(s')\right\} Q\mathcal{L}(s) P\mathcal{W}(s). \tag{19}$$

Inserting Equation (19) in Equation (17) and using the properties of *P* we ge<sup>t</sup> the Nakajima–Zwanzig equation:

$$\begin{split} i\hbar \boldsymbol{P} \frac{\partial \mathcal{W}(t)}{\partial t} &= \boldsymbol{P} \mathcal{L}\_{S} \boldsymbol{W}(t) \\ &+ \quad \frac{1}{i\hbar} \boldsymbol{P} \mathcal{L}\_{SR}(t) \boldsymbol{Q} \int\_{t\_{0}}^{t} ds \boldsymbol{T} \exp\left\{-\frac{i}{\hbar} \int\_{s}^{t} ds' \boldsymbol{Q} \mathcal{L}(s') \boldsymbol{Q}\right\} \boldsymbol{Q} \mathcal{L}\_{SR}(s) \boldsymbol{P} \mathcal{W}(s). \end{split} \tag{20}$$

In order to have an explicit perturbative expansion in powers of *HSR*(*t*) one has to factorize the time-ordered exponential as follows:

$$T \exp\left\{-\frac{i}{\hbar} \int\_s^t ds' Q \mathcal{L}(s') Q\right\} = \exp\{Q \mathcal{L}\_0 Q\} (1 + \mathcal{R}),\tag{21}$$

where the remainder 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 *HSR*:

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

Now, by taking into account that for any operator *A* acting on the Fock space of the hybrid system *e* −*it*L<sup>0</sup> *A* = *e*<sup>−</sup>*itH*<sup>0</sup> *AeitH*0 and denoting by *<sup>U</sup>*0(*<sup>t</sup>*,*<sup>s</sup>*) = *e*<sup>−</sup>*<sup>i</sup>*(*<sup>t</sup>*−*<sup>s</sup>*)*<sup>H</sup>*<sup>0</sup> the unitary evolution of the disconnected systems we arrive at the well known form of the GME:

$$\begin{split} \dot{\boldsymbol{\rho}}(t) &= \begin{bmatrix} \boldsymbol{H}\_{\mathcal{S}\prime}\boldsymbol{\rho}(t) \end{bmatrix} - \frac{i}{\hbar} \boldsymbol{l} \boldsymbol{l}\_{0}^{\dagger}(t, t\_{0}) \text{Tr}\_{\mathcal{R}} \left\{ \int\_{t\_{0}}^{t} \boldsymbol{ds} \left[ \boldsymbol{H}\_{\mathcal{S}\prime}(t)\_{\prime} \left[ \boldsymbol{H}\_{\mathcal{S}\prime}(\boldsymbol{s})\_{\prime} \boldsymbol{\tilde{\rho}}(\boldsymbol{s}) \boldsymbol{\rho}\_{\mathcal{R}} \right] \right] \right\} \boldsymbol{l} \boldsymbol{l}\_{0}(t, t\_{0}) \\ &= \begin{bmatrix} \boldsymbol{H}\_{\mathcal{S}\prime}\boldsymbol{\rho}(t) \end{bmatrix} - \frac{i}{\hbar} \boldsymbol{l} \boldsymbol{l}\_{\mathcal{S}}^{\dagger}(t, t\_{0}) \text{Tr}\_{\mathcal{R}} \left\{ \int\_{t\_{0}}^{t} \boldsymbol{ds} \left[ \boldsymbol{\tilde{\rho}}\_{\mathcal{S}\prime}(t)\_{\prime} \left[ \boldsymbol{\tilde{\rho}}\_{\mathcal{S}\prime}(\boldsymbol{s})\_{\prime} \boldsymbol{\tilde{\rho}}(\boldsymbol{s}) \boldsymbol{\rho}\_{\mathcal{R}} \right] \right] \right\} \boldsymbol{l} \boldsymbol{l}\_{\mathcal{S}}(t, t\_{0}), \end{split}$$

where in order to ge<sup>t</sup> 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 *HT* and *HE* carries only one creation or annihilation operator for the corresponding reservoir such that:

$$\text{Tr}\_{\mathbb{R}}\{\tilde{c}\_{ql}^{\dagger}(t)\tilde{b}\_{\mathbf{k}}(s)\rho\_{\mathcal{R}}\} = \text{Tr}\_{\text{leading}}\{\tilde{c}\_{ql}^{\dagger}(t)\rho\_{\text{levels}}\} \cdot \text{Tr}\_{\text{bath}}\{\tilde{b}\_{\mathbf{k}}(s)\rho\_{\text{bath}}\} = 0. \tag{23}$$

Moreover, the time evolution of each term can be simplified due to the commutation relations [*<sup>H</sup>*bath, *HT*]=[*<sup>H</sup>*leads, *HE*] = 0:

$$H\_T(t) \quad = \; e^{\dot{\xi}tH\_S} e^{\dot{\xi}tH\_{\text{laus}}} H\_T e^{-\dot{\xi}tH\_S} e^{-\dot{\xi}tH\_{\text{laus}}} \; \tag{24}$$

$$
\bar{H}\_{\rm E}(t) \quad = \; e^{\dot{\xi}tH\_{\rm S}}e^{\dot{\xi}tH\_{\rm bath}}H\_{\rm E}e^{-\dot{\xi}tH\_{\rm S}}e^{-\dot{\xi}tH\_{\rm bath}}.\tag{25}
$$

The GME then reads as:

$$\begin{split} \rho(t) &= -\frac{i}{\hbar}[H\_S,\rho(t)] - \frac{1}{\hbar^2} \mathcal{U}\_S^\dagger(t, t\_0) \text{Tr}\_{\text{leading}} \left\{ \int\_{t\_0}^t ds \left[ \bar{H}\_T(t), \left[ \bar{H}\_T(s), \bar{\rho}(s) \rho\_{\text{back}} \right] \right] \right\} \mathcal{U}\_S(t, t\_0) \\ &- \frac{1}{\hbar^2} \mathcal{U}\_S^\dagger(t, t\_0) \text{Tr}\_{\text{bath}} \left\{ \int\_{t\_0}^t ds \left[ \bar{H}\_E(t), \left[ \bar{H}\_E(s), \bar{\rho}(s) \rho\_{\text{back}} \right] \right] \right\} \mathcal{U}\_S(t, t\_0) \end{split} \tag{26}$$

$$\dot{\rho} := -\frac{i}{\hbar} [H\_{S\prime} \rho(t)] - \mathcal{D}\_{\text{hads}}[\rho\_{\prime} t] - \mathcal{D}\_{\text{bath}}[\rho\_{\prime} t]. \tag{27}$$

Equation (27) is the generalized master equation for our hybrid system. It provides the reduced density operator *ρ* 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 *χl*. The third term in Equation (27) is needed only if *HS*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 *gEM* (see Section 5). On the other hand, for our calculations *gEM h*¯ *ω*, *ω* being the frequency of the cavity mode. Then the RWA holds and <sup>D</sup>bath[*ρ*, *t*] can be cast in a Lindblad form. Let us stress that in the ultrastrong coupling regime on typically has *gEM*/¯*hω* > 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 F*S* = F*S*1⊗ F*S*2. We discuss this issue in the next subsection.
