*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 *HT* and *HE* w.r.t. the 'disjointed' basis formed by the eigenstates of *HS*1 and *HS*2 , that is |*ν*, *j* := |*ν*⊗|*j*. However this strategy does not help much when evaluating the time evolution (see Equations (24) and (25)) as *HS* is not diagonal w.r.t. to |*ν*, *j* such that one cannot easily write down the matrix elements of the unitary evolution *US*(*<sup>t</sup>*, *<sup>t</sup>*0). In fact we are forced to solve the GME by using the eigenstates |*ϕp*) and eigenvalues E*p* of the Hamiltonian *HS*. The former are written as:

$$\mathbb{V}\left(\boldsymbol{\varphi}\_{\mathcal{P}}\right) = \sum\_{\boldsymbol{\nu}\_{\mathcal{I}}\mathbf{j}} \mathcal{V}^{(p)}\_{\boldsymbol{\nu}\mathbf{j}} \left| \boldsymbol{\nu}\_{\boldsymbol{\nu}} \boldsymbol{j} \right\rangle. \tag{28}$$

Here the round bracket notation |*ϕp*) is meant to underline that the state *ϕ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 *HS*. This notation also prevents any confusion if the 'free' states |*ν*, *j* were also labelled by a single index *p*. In that case the above equation is conveniently rewritten as |*ϕp*) = ∑*p* V(*p*) *p* |*p*. 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 |*ν*, *j* and is not necessarily small. Therefore we shall not follow a perturbative approach but rather calculate E*p* and the weights V(*p*) *νj* by numerically diagonalizing *HS* 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 *US* describes the joint systems *S*1 and *S*2 and therefore the hybrid interaction cannot be simply neglected neither in Dleads nor in Dbath; in fact one can easily check that *V* does not commute with *HE* or *HT*. Moreover, as has been clearly pointed out by Beaudoin et al. [30], by disregarding the qubit-resonator interaction when calculating Dbath 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 *HS* poses serious technical problems because both spaces F*S*1 and F*S*2 are in principle infinite dimensional. Besides that, the Coulomb interaction in *HS*1 prevents one to derive the interacting many-body configurations {|*ν*} 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 ˆ *h*(0) *S*1 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*ses single-particle states {*ψ*1, *ψ*2, ..., *ψ <sup>N</sup>*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*mes non-interacting many-body configurations (NMBS) {| *<sup>λ</sup>*}*<sup>λ</sup>*=1,..., *N*mes from the *N*ses single-particle states introduced above. Note that for computational reasons we have to keep *N*mes < 2*N*ses for larger *N*ses. Then, if *i λ n* is the occupation number of the single-particle state *ψ<sup>n</sup>*, a non-interacting many-body configuration |*λ* reads as:

$$|\lambda\rangle = |i\_1^{\lambda}, i\_1^{\lambda}, \dots, i\_{N\_{\text{sum}}}^{\lambda}\rangle. \tag{29}$$

(*D*3) Diagonalization of the Coulomb-interacting electronic Hamiltonian *HS*1 = *H*(0) *S*1 + *W* on the subspace of non-interacting many-body states from F*S*1 . As a result one gets *N*mes interacting many-body states (IMBS) |*ν* and the associated energy levels *Eν* introduced in Section 2.1. We also introduce the 'free' energies of *HS*1 + *HS*2 , that is E(0) *<sup>ν</sup>*,*j* = *Eν* + *ej*. Note that in view of diagonalization the interaction *V* between the two subsystems must be also written w.r.t. the 'free' states {|*<sup>ν</sup>*, *j*}. If the second subsystem is not needed then the GME must be solved w.r.t. the set {|*ν*} 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 *US* is diagonal w.r.t. the many-body basis {|*ν*}, the matrix elements of the fermionic operators in the interaction picture *ν*|*c*˜*nσ*(*t*)|*ν* depend on energy differences *Eν* − *Eν* which, due to the Coulomb effects, cannot be reduced to the single-particle energy *εn*.

(*D*4) Diagonalization of the fully interacting Hamiltonian *HS* on a subspace of F made by the lowest energy *N*mesT interacting MBS of *HS*1 and *j*max eigenstates of *HS*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*mesT < *N*mes) and w.r.t. the states in F*S*2.

Once this procedure is performed, one can express the system-environment couplings *HT* and *HB* in the fully interacting basis and use the eigenvalues E*p* to replace the unitary evolution *US* by the corresponding diagonal matrix *e*<sup>−</sup>*it*E*<sup>p</sup> <sup>δ</sup>pp* . 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 *HS*.

(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*mes*<sup>T</sup>* interacting MBSs are derived from a larger set of non-interacting states {| *<sup>λ</sup>*}*<sup>λ</sup>*=1,..., *N*mes . Then the calculated *Np* := *N*mes*<sup>T</sup>* × *j*max fully interacting states are more accurate than the ones obtained by diagonalizing once a *Np* × *Np* matrix. On the other hand, enlarging the full space to *N*mes × *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*ses single-electron states we construct the set of non-interacting MBSs containing up to *Ne* electrons, the size of this set being, of course, ( *N*ses *Ne* ). The accuracy of numerical diagonalization which leads to the interacting many-body configurations with up to *Ne* electrons is essentially assessed by comparing the spectra associated obtained for different *N*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 ge<sup>t</sup> stable transport simulations.

(b) Physical interpretation. It is obvious that the Coulomb interaction *W* mixes only the non-interacting many-body configurations |*λ* while the hybrid coupling *V* mixes both *λ* and *j* states. For this reason the weights of a non-interacting state |*λ*, *j* in a fully interacting state |*ϕ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 *HS*. 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 ye<sup>t</sup> 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<sup>10</sup> 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.
