*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\_a - a)} \approx \pi \delta(E\_\beta - E\_a - \epsilon),\tag{60}$$

where a small imaginary principle part is ignored. We have furthermore assumed instant lead-system coupling at *t* = 0 with *<sup>χ</sup>l*(*t*) = *<sup>θ</sup>*(*t*), the Heaviside unit step function, in Equation (8) for *HT*. 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

$$\text{vec}(A\rho B) = (B^T \otimes A)\text{vec}(\rho) \tag{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* ⊗ *A* = {*<sup>B</sup><sup>α</sup>*,*β<sup>A</sup>*}. The Kronecker product of two *N*mes × *N*mes matrices results in a *N*2mes × *N*2mes 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–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\_l \text{vec}(\rho) = \mathcal{L} \text{vec}(\rho),
\tag{62}
$$

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

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

in terms of the matrices of the column stacked left U, and the right V eigenvectors of L

$$
\mathcal{L}\mathcal{V} = \mathcal{V}\mathcal{L}\_{\text{diag}} \quad \text{and} \quad \mathcal{U}\mathcal{L} = \mathcal{L}\_{\text{diag}}\mathcal{U}.\tag{64}
$$

obeying

$$
\mathcal{U}\mathcal{V} = I,\quad\text{and}\quad\mathcal{V}\mathcal{U} = I.\tag{65}
$$

Special care has to be taken in the numerical implementation of this solution procedure as many software packages use another normalization for U and 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 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(\tau) = \text{Tr}\_{\mathbb{R}}\left\{ e^{-iH\tau/\hbar} \mathcal{X} \rho\_{\mathbb{T}}(0) e^{+iH\tau/\hbar} \right\},\tag{66}$$

where *H* is the total Hamiltonian of the system, *ρ*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*Λ*<sup>l</sup>* for current correlations as in [107], where *Q* = ∑*i c*†*i ci* is the fermionic charge operator and Λ*<sup>l</sup>* is the Liouville dissipation operator for lead *l*. The structure of *χ* (66) indicates that the two-time correlation function is then

$$
\langle X(\tau)X(0) \rangle = \text{Tr}\_{\mathbb{S}} \left\{ X(0) \chi(\tau) \right\},
\tag{67}
$$

with

$$\text{vec}(\chi(\tau)) = \{\mathcal{U}[\exp\left(\mathcal{L}\_{\text{diag}}t\right)]\mathcal{V}\}\text{vec}(\chi(0)).\tag{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 *<sup>S</sup>*(*E*), and for the current-current correlation the corresponding Fourier spectral density denoted by *Dll*(*E*), where *l* and *l* refer to *L* and *R*, the Left and Right leads.
