Markovian treatment of non-Markovian dynamics of open Fermionic systems

Open nonequilibirum systems are at the forefront of experimental and theoretical research due to the rich and complex physics they provide access to as well as due to applicational prospects of building nanoscale devices for quantum based technologies and computations [1–3]. Especially intriguing in term of both fundamental science and potential applications are effects of strong correlations. A number of impurity solvers capable of treating strongly correlated systems coupled to continuum of baths degrees of freedom were developed. Among them are numerical renormalization group in the basis of scattering states [4, 5], flow equations [6, 7], time-dependent density matrix renormalization group [8, 9], multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) [10, 11], and continuous time quantum Monte Carlo [12–14] approaches. These numerically exact techniques are very demanding and so far are mostly applicable to simple models only.

At the same time, accurate numerically inexpensive impurity solvers are in great demand both as standalone techniques to be applied in simulation of, e.g. nanoscale junctions and as a part of divide-and-conquer schemes such as, e.g. dynamical mean-field theory (DMFT) [15, 16]. In this respect ability to map complicated non-Markovian dynamics of a system onto much simpler Markov consideration is an important step towards creating new computational techniques applicable in realistic simulations. In particular, such mapping was used in auxiliary master equation approach (AMEA) [17, 18] introducing numerically inexpensive and pretty accurate solver for the nonequilibrium DMFT. Within AMEA the original unitary evolution is substituted with a Lindblad-type quantum master equation consideration of an expanded system (system plus set of auxiliary modes). Another example is the recent formulation of the auxiliary dual-fermion (aux-DF) method [19]. Aux-DF technique complements the mapping with a procedure correcting for deviation of original hybridization functions from those in the auxiliary system. We note that similar mapping ideas are employed also in the reaction coordinate formalism [20–22]. While the mappings appear to be very useful and accurate, in most cases only semi-quantitative arguments to justify the mapping were presented with main supporting evidence being benchmarking versus numerically exact computational techniques. In particular, a justification for the mapping was put forward in [23–25] based upon the singular coupling derivation of the Lindblad equation.

Recently, a rigorous proof of non-Markov to Markov mapping for open Bose quantum systems was presented in the literature [26]. It was shown that the evolution of reduced density matrix in non-Markov system with unitary system-environment evolution can be equivalently obtained by a Markov evolution of an extended system (system plus modes of environment) under non-unitary (Lindblad-type) evolution. Here, we extend the consideration of [26] to fermionic open quantum systems and to multi-time correlation functions. The structure of the paper is the following. After introducing physical and auxiliary models of an open quantum Fermi system in section 2 we discuss non-Markov to Markov mapping procedure in section 3. Exact mathematical proofs are given in Appendices. Section 4 presents numerical illustration of the mapping for a simple generic model of a junction. We conclude in section 5.

We consider an open fermionic system S coupled to an arbitrary number N of external baths, initially each at its own thermodynamic equilibrium, i.e. characterized by its own electrochemical potential and temperature (see figure 1(a)). The Hamiltonian of the model is

$$\begin{eqnarray}&&{\hat{H}}^{{\rm{phys}}}(t)={\hat{H}}_{S}(t)+\displaystyle \sum _{B=1}^{N}\left({\hat{H}}_{B}+{\hat{V}}_{{SB}}\right).\end{eqnarray} \tag{ 1 }$$

Here ${\hat{H}}_{S}(t)$ and ${\hat{H}}_{B}$ (B ∈ {1, ..., N}) are Hamiltonians of the system and baths. ${\hat{V}}_{{SB}}$ introduces coupling of the system to bath B. While the Hamiltonian of the system is general and may be time-dependent, we follow the usual paradigm by assuming bi-linear coupling in constructing fermionic junction models

$$\begin{eqnarray}&&{\hat{H}}_{B}=\displaystyle \sum _{k\in B}{\varepsilon }_{{Bk}}\,{\hat{c}}_{{Bk}}^{\dagger }{\hat{c}}_{{Bk}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\hat{V}}_{{SB}}=\displaystyle \sum _{k\in B}\displaystyle \sum _{i\in S}({V}_{i,{Bk}}\,{\hat{d}}_{i}^{\dagger }{\hat{c}}_{{Bk}}+{\rm{h}}.{\rm{c}}.),\end{eqnarray} \tag{ 3 }$$

where ${\hat{d}}_{i}^{\dagger }$ (${\hat{d}}_{i}$) and ${\hat{c}}_{{Bk}}^{\dagger }$ (${\hat{c}}_{{Bk}}$) create (annihilate) electron in level i of the system S and level k of bath B. In the model, dynamics of the system-plus-baths evolution is unitary. Below we call this model phys (physical). We note in passing that extension of the consideration to other types of system–baths couplings is straightforward, as long as baths are quadratic in the Fermi operators.

Zoom In Zoom Out Reset image size

The other configuration we will consider is a model we shall call aux (auxiliary; see figure 1(b)). Here, the same system S is coupled to a number of auxiliary modes A, which in their turn are coupled to two baths. There are two Fermi baths in the configuration: one (L) is completely full (${\mu }_{L}\to +\infty$), the other (R) is completely empty (${\mu }_{R}\to -\infty$). The Hamiltonian of the total system is

$$\begin{eqnarray}&&{\hat{H}}^{{\rm{aux}}}(t)={\hat{H}}_{S}(t)+{\hat{V}}_{{SA}}+{\hat{H}}_{A}+\displaystyle \sum _{C=L,R}\left({\hat{H}}_{C}+{\hat{V}}_{{AC}}\right),\end{eqnarray} \tag{ 4 }$$

where ${\hat{H}}_{S}$ is the same as in (1), ${\hat{H}}_{A}$ represents set of modes

$$\begin{eqnarray}&&{\hat{H}}_{A}=\displaystyle \sum _{{m}_{1},{m}_{2}\in A}{H}_{{m}_{1}{m}_{2}}^{A}{\hat{a}}_{{m}_{1}}^{\dagger }{\hat{a}}_{{m}_{2}}\end{eqnarray} \tag{ 5 }$$

and ${\hat{V}}_{{SA}}$ their interaction with the system

$$\begin{eqnarray}&&{\hat{V}}_{{SA}}=\displaystyle \sum _{i\in S}\displaystyle \sum _{m\in A}({V}_{{im}}^{{SA}}{\hat{d}}_{i}^{\dagger }{\hat{a}}_{m}+{\rm{h}}{\rm{.}}{\rm{c}}.).\end{eqnarray} \tag{ 6 }$$

Here ${\hat{a}}_{m}^{\dagger }$ (${\hat{a}}_{m}$) creates (annihilates) electron in the auxiliary mode m in A.

${\hat{H}}_{C}$ represents continuum of states in contact C

$$\begin{eqnarray}&&{\hat{H}}_{C}=\displaystyle \sum _{k\in C}{\varepsilon }_{{Ck}}{\hat{c}}_{{Ck}}^{\dagger }{\hat{c}}_{{Ck}}\end{eqnarray} \tag{ 7 }$$

with constant density of states

$$\begin{eqnarray}&&{N}_{C}(E)\equiv \displaystyle \sum _{k\in C}\delta (E-{\varepsilon }_{{Ck}})={\rm{const}}\end{eqnarray} \tag{ 8 }$$

and ${\hat{V}}_{{AC}}$ couples auxiliary modes A to bath C (L or R)

$$\begin{eqnarray}&&{\hat{V}}_{{AC}}=\displaystyle \sum _{k\in C}\displaystyle \sum _{m\in A}({t}_{m}^{C}{\hat{a}}_{m}^{\dagger }{\hat{c}}_{{Ck}}+{\rm{h}}{\rm{.}}{\rm{c}}.).\end{eqnarray} \tag{ 9 }$$

Dynamics of the whole configuration is unitary.

In the next section we show that the reduced time evolution of S in models phys and aux is the same (subject to certain conditions) and that the reduced dynamics of S + A in model aux satisfies an appropriate Lindblad Markov evolution. This establishes procedure for Markov non-unitary Lindblad-type treatment of S + A in aux exactly representing overall (i.e. system plus baths) unitary non-Markov dynamics of S in phys by tracing out A degrees of freedom.

Consideration of the mapping consists of three levels of description: (1) overall (S plus baths) unitary dynamics of the physical system (phys); (2) overall (S + A plus baths) unitary dynamics of the auxiliary system (aux); (3) non-unitary Lindblad-type dynamics of S + A in the auxiliary system (aux). Below we first discuss equivalence of the unitary dynamics of S in phys and aux systems, then we prove equivalence of unitary and Lindblad-type evolutions in the aux system.

First, we are going to prove that with an appropriate choice of parameters of aux the dynamics of S can be equivalently represented in the original model phys and auxiliary model aux, under assumption that the dynamics of the whole system is unitary. Because non-interacting baths are fully characterized by their two-time correlation functions, equivalence of system-bath(s) hybridizations (i.e. correlation functions of the bath(s) dressed with system-bath(s) interactions) for the two models indicates equivalence of the reduced system dynamics in the two cases. For example, hybridization function is the only information about baths in numerically exact simulations of strongly correlated systems [13]. Nonequilibrium character of the system requires fitting two projections of the hybridization function (also called self-energy in the literature). In particular, these may be retarded and Keldysh projections. Let ${{\rm{\Sigma }}}_{B}^{r}(E)$ and ${{\rm{\Sigma }}}_{B}^{K}(E)$ be matrices introducing the corresponding hybridization functions for bath B of the physical problem (figure 1(a)). Retarded projection carries information on bath's spectral function and strength of system-bath coupling yielding dissipation rates for the system due to coupling to the bath. Keldysh projection yields information on bath's population which by Pauli principle defines possibility of electron exchange between system and bath

$$\begin{eqnarray}&&{\left({{\rm{\Sigma }}}_{B}^{r}(E)\right)}_{{ij}}=\displaystyle \sum _{k\in B}{V}_{i,{Bk}}\,{g}_{{Bk}}^{r}(E){V}_{{Bk},j},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\left({{\rm{\Sigma }}}_{B}^{K}(E)\right)}_{{ij}}=\displaystyle \sum _{k\in B}{V}_{i,{Bk}}\,{g}_{{Bk}}^{K}(E){V}_{{Bk},j},\end{eqnarray} \tag{ 11 }$$

where ${g}_{{Bk}}^{r(K)}(E)$ are the Fourier transforms of retarded (Keldysh) projections of the decoupled (${\hat{V}}_{{SB}}=0$) electron Green'(s) function ${g}_{{Bk}}(\tau ,\tau ^{\prime} )=-{\rm{i}}\langle {T}_{c}\,{\hat{c}}_{{Bk}}(\tau ){\hat{c}}_{{Bk}}^{\dagger }(\tau ^{\prime} )\rangle$ in contact B. The parameters of the auxiliary model should then be chosen such that the total hybridization functions for the system

$$\begin{eqnarray}&&{{\rm{\Sigma }}}^{r}(E)=\displaystyle \sum _{B=1}^{N}{{\rm{\Sigma }}}_{B}^{r}(E),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}^{K}(E)=2\,{\rm{i}}\displaystyle \sum _{B=1}^{N}\left(1-2{f}_{B}(E)\right){\rm{Im}}\,{{\rm{\Sigma }}}_{B}^{r}(E)\end{eqnarray} \tag{ 13 }$$

are as close as possible to the corresponding hybridization functions, ${\tilde{{\rm{\Sigma }}}}^{r}(E)$ and ${\tilde{{\rm{\Sigma }}}}^{K}(E)$, of S in the auxiliary model (figure 1(b)) [17, 18, 24]. The latter have contribution from full (L) and empty (R) baths, and from auxiliary modes (A)

$$\begin{eqnarray}&&{\tilde{{\rm{\Sigma }}}}^{r}(E)={\tilde{{\rm{\Sigma }}}}_{L}^{r}(E)+{\tilde{{\rm{\Sigma }}}}_{R}^{r}(E),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\tilde{{\rm{\Sigma }}}}^{K}(E)=2\,i\,{\rm{Im}}\left({\tilde{{\rm{\Sigma }}}}_{R}^{r}(E)-{\tilde{{\rm{\Sigma }}}}_{L}^{r}(E)\right),\end{eqnarray} \tag{ 15 }$$

where we assume modes A initially in equilibrium with its contact (or in stationary state if coupled to L and R). Requirement of equivalence can be expressed as

$$\begin{eqnarray}&&{\rm{Im}}\,{\tilde{{\rm{\Sigma }}}}_{L}^{r}(E)=\displaystyle \frac{2\,{\rm{i}}\,{\rm{Im}}{{\rm{\Sigma }}}^{r}(E)+{{\rm{\Sigma }}}^{K}(E)}{4\,i},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\rm{Im}}\,{\tilde{{\rm{\Sigma }}}}_{R}^{r}(E)=\displaystyle \frac{2\,{\rm{i}}\,{\rm{Im}}{{\rm{\Sigma }}}^{r}(E)-{{\rm{\Sigma }}}^{K}(E)}{4\,i}.\end{eqnarray} \tag{ 17 }$$

Thus, the problem reduces to fitting known functions in the right side of the expression with multiple contributions from auxiliary modes to the hybridization functions in the left side. We note that the knowledge of total (sum of contributions from all baths) hybridization function (retarded and Keldysh components) allows to fully determine interacting correlation functions in the S subspace of phys. That is, no information on contribution from each separate bath is required. Thus, any number of baths B in physical system can be represented by only two baths (one full and one empty) in the auxiliary system. The exact mapping we prove below allows to evaluate correlation functions (and in particular, single particle Green's functions) in the S subspace of the physical system by considering Lindblad-type evolution in the aux system. After the correlation functions has been evaluated fluxes between the system S and baths B can be evaluated utilizing the well-known exact Jauho–Meir–Wingreen and similar expressions.

In principle fitting (16), (17) can be done in many different ways [24]. For example, possibility of exact fitting of an arbitrary function with set of Lorentzians was discussed in [27]. In auxiliary systems such fitting corresponds to a construction where each auxiliary mode is coupled to its own bath. Note that in practical simulations accuracy of the results can be improved either by increasing number of auxiliary modes or by considering more general (nondiagonal) level-bath geometries in the auxiliary system, as is implemented in, e.g. AMEA [28], or by employing diagrammatic expansion related to the difference between true and fitted hybridization functions, as is realized in, e.g. dual fermion approach [29], or both.

After having established the equivalence of reduced system (S) dynamics in phys and aux, we now turn to consideration of evolution of the aux model. We will show that reduced S + A dynamics derived from unitary evolution of the aux model can be exactly represented by a suitable non-unitary Lindblad-type evolution.

Following [26] we consider the reduced density operator of S + A in aux, ${\hat{\rho }}_{{SA}}$, which is defined by integrating out the baths degrees of freedom from the total density operator ${\hat{\rho }}^{{\rm{aux}}}(t)$

$$\begin{eqnarray}&&{\hat{\rho }}_{{SA}}(t)\equiv {{\rm{Tr}}}_{{LR}}\left[{\hat{\rho }}^{{\rm{aux}}}(t)\right]\end{eqnarray} \tag{ 18 }$$

${\hat{\rho }}^{{\rm{aux}}}$ follows an unitary evolution with initial condition given by S + A decoupled from the baths

$$\begin{eqnarray}&&{\hat{\rho }}^{{\rm{aux}}}(0)={\hat{\rho }}_{L}\otimes {\hat{\rho }}_{{SA}}(0)\otimes {\hat{\rho }}_{R},\end{eqnarray} \tag{ 19 }$$

where ${\hat{\rho }}_{L}=| {\rm{full}}\rangle \langle {\rm{full}}|$, ${\hat{\rho }}_{R}=| {\rm{empty}}\rangle \langle {\rm{empty}}|$, and ${\hat{\rho }}_{{SA}}(0)$ arbitrary.

In appendix A we prove that ${\hat{\rho }}_{{SA}}(t)$ satisfies the following Markov Lindblad-type equation of motion

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\hat{\rho }}_{{SA}}(t) & = & -{\rm{i}}\left[{\hat{H}}_{{SA}}(t),{\hat{\rho }}_{{SA}}(t)\right]\\ & & +\ \displaystyle \sum _{{m}_{1},{m}_{2}\in A}\left[{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}\left({\hat{a}}_{{m}_{2}}{\hat{\rho }}_{{SA}}(t){\hat{a}}_{{m}_{1}}^{\dagger }-\displaystyle \frac{1}{2}\left\{{\hat{\rho }}_{{SA}}(t),{\hat{a}}_{{m}_{1}}^{\dagger }{\hat{a}}_{{m}_{2}}\right\}\right)\right.\\ & & +\ \left.{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}\left({\hat{a}}_{{m}_{1}}^{\dagger }{\hat{\rho }}_{{SA}}(t){\hat{a}}_{{m}_{2}}-\displaystyle \frac{1}{2}\left\{{\hat{\rho }}_{{SA}}(t),{\hat{a}}_{{m}_{2}}{\hat{a}}_{{m}_{1}}^{\dagger }\right\}\right)\right]\\ & \equiv & \ {{ \mathcal L }}_{{SA}}(t)| {\rho }_{{SA}}(t)\rangle \rangle ,\end{array}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{\hat{H}}_{{SA}}(t)\equiv {\hat{H}}_{S}(t)+{\hat{V}}_{{SA}}+{\hat{H}}_{A},\end{eqnarray} \tag{ 21 }$$

${{ \mathcal L }}_{{SA}}$ is the Liouvillian superoperator defined on the S + A subspace of the aux model and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{C}\equiv 2\pi {t}_{{m}_{1}}^{C}{\left({t}_{{m}_{2}}^{C}\right)}^{* }{N}_{C}\ (C=L,R)\end{eqnarray} \tag{ 22 }$$

is the dissipation matrix due to the coupling to contact C.

Next we turn to multi-time correlation functions of operators in the S + A subspace of the aux model. Following [26] we start consideration from two-time correlation function on real time axis. For arbitrary operators ${\hat{O}}_{1}$ and ${\hat{O}}_{2}$ in S + A we define two-time (t₁ ≥ t₂ ≥ 0) correlation function as

$$\begin{eqnarray}&&\langle {\hat{O}}_{1}({t}_{1}){\hat{O}}_{2}({t}_{2})\rangle \equiv {\rm{Tr}}\left[{\hat{O}}_{1}\,{\hat{U}}^{{\rm{aux}}}({t}_{1},{t}_{2}){\hat{O}}_{2}\,{\hat{U}}^{{\rm{aux}}}({t}_{2},0){\hat{\rho }}^{{\rm{aux}}}(0){\hat{U}}^{{\rm{aux}}\dagger }({t}_{1},0)\right].\end{eqnarray} \tag{ 23 }$$

Here ${\hat{U}}^{{\rm{aux}}}$ is the evolution operator in the aux system

$$\begin{eqnarray}&&{\hat{U}}^{{\rm{aux}}}(t,t^{\prime} )\equiv T\exp \left[-{\rm{i}}{\int }_{t^{\prime} }^{t}{\rm{d}}{s}\,{\hat{H}}^{{\rm{aux}}}(s)\right]\end{eqnarray} \tag{ 24 }$$

and T is the time-ordering operator. In appendix B we show that (23) can be equivalently obtained from reduced Lindblad-type evolution in the S + A subspace

$$\begin{eqnarray}&&\langle {\hat{O}}_{1}({t}_{1}){\hat{O}}_{2}({t}_{2})\rangle =\langle \langle I| {{ \mathcal O }}_{1}^{-}\,{{ \mathcal U }}_{{SA}}({t}_{1},{t}_{2}){{ \mathcal O }}_{2}^{-}\,{{ \mathcal U }}_{{SA}}({t}_{2},0)| {\rho }_{{SA}}(0)\rangle \rangle .\end{eqnarray} \tag{ 25 }$$

Here $\langle \langle I|$ is Liouville space bra representation of the Hilbert space identity operator, $| {\rho }_{{SA}}(0)\rangle \rangle$ is Liouville space ket representation of the Hilbert space operator ${\hat{\rho }}_{{SA}}(0)$, ${{ \mathcal O }}_{i}$ is the Liouville space superoperator corresponding to the Hilbert space operator ${\hat{O}}_{i}$ (see figure 2)

$$\begin{eqnarray}{{ \mathcal O }}_{i}| \rho \rangle \rangle =\left\{\begin{array}{cc}{{ \mathcal O }}_{i}^{-}| \rho \rangle \rangle \equiv {\hat{O}}_{i}\,\hat{\rho } & {\rm{forward}}\,{\rm{branch}}\\ {{ \mathcal O }}_{i}^{+}| \rho \rangle \rangle \equiv \hat{\rho }\,{\hat{O}}_{i} & {\rm{backward}}\,{\rm{branch}}\end{array}\right.\end{eqnarray} \tag{ 26 }$$

and ${{ \mathcal U }}_{{SA}}$ is the Liouville space evolution superoperator

$$\begin{eqnarray}&&{{ \mathcal U }}_{{SA}}(t,t^{\prime} )\equiv T\exp \left[{\int }_{t^{\prime} }^{t}{\rm{d}}{s}\,{{ \mathcal L }}_{{SA}}(s)\right].\end{eqnarray} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

Finally, we extend consideration to multi-time correlation functions of arbitrary operators ${\hat{O}}_{i}$ (i ∈ {1, ..., N}) defined on the Keldysh contour (see figure 2) as

$$\begin{eqnarray}&&\langle {T}_{c}\,{\hat{O}}_{1}({\tau }_{1}){\hat{O}}_{2}({\tau }_{2})\ldots {\hat{O}}_{N}({\tau }_{N})\rangle \equiv {\rm{Tr}}\left[{T}_{c}\,{\hat{O}}_{1}\,{\hat{O}}_{2}\ldots {\hat{O}}_{N}\,{\hat{U}}_{c}\,{\hat{\rho }}^{{\rm{aux}}}(0)\right],\end{eqnarray} \tag{ 28 }$$

where τ_i are the contour variables, T_c is the contour ordering operator, and

$$\begin{eqnarray}&&{\hat{U}}_{c}={T}_{c}\exp \left[-{\rm{i}}{\int }_{c}{\rm{d}}\tau \,{\hat{H}}^{{\rm{aux}}}(\tau )\right]\end{eqnarray} \tag{ 29 }$$

is the contour evolution operator. Note subscripts of operators O_i in the right side of (28) indicate operators on either side of the contour. In appendix C we prove that multi-time correlation functions (28) can be evaluated solely from Markov Lindblad-type evolution in S + A subspace of the aux model

$$\begin{eqnarray}&&\begin{array}{l}\langle {T}_{c}\,{\hat{O}}_{1}({\tau }_{1}){\hat{O}}_{2}({\tau }_{2})\ldots {\hat{O}}_{N}({\tau }_{N})\rangle \\ \ =\ {\left(-1\right)}^{P}\langle \langle I| {{ \mathcal O }}_{{\theta }_{1}}\,{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{1}},{t}_{{\theta }_{2}}){{ \mathcal O }}_{{\theta }_{2}}\,{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{2}},{t}_{{\theta }_{3}})\ldots {{ \mathcal O }}_{{\theta }_{N}}{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{N}},0)| {\rho }_{{SA}}(0)\rangle \rangle .\end{array}\end{eqnarray} \tag{ 30 }$$

Here P is number of Fermi interchanges in the permutation of operators ${\hat{O}}_{i}$ by T_c, θ_i are indices of operators ${\hat{O}}_{i}$ rearranged is such a way that${t}_{{\theta }_{1}}\gt {t}_{{\theta }_{2}}\gt \mathrm{...}\gt {t}_{{\theta }_{N}}$ (${t}_{{\theta }_{i}}$ is real time corresponding to contour variable ${\tau }_{{\theta }_{i}}$), ${{ \mathcal O }}_{{\theta }_{i}}$ are the superoperators defined in (26), and ${{ \mathcal U }}_{{SA}}$ is the Liouville space evolution superoperator defined in (27).

Equivalence of S dynamics derived from unitary evolution of models phys and aux together with (20) and (30) completes the proof of the possibility of Markov treatment for non-Markovian dynamics in open quantum Fermi systems.

Application of the mapping in realistic simulations relies on ability to fit hybridization function of the phys system with a set of auxiliary modes in the aux system, equations (16), (17). In general, to fit arbitrary function one needs infinite number of auxiliary modes, while in realistic calculations one can account for only final number of modes. Thus, when applying the mapping one is looking for a trade-off between accuracy and efficiency: the more auxiliary modes are considered the better is the fit and the more involved is procedure to solve the auxiliary quantum master equations (QME). For example, in [28] high accuracy of fitting with 16 auxiliary modes was demonstrated for the Anderson impurity model. However, for aux system of this size application of the matrix product states was necessary to solve the QME.

Alternative (or complementary) strategy is utilization of small number of auxiliary modes (thus, relying on relatively poor fitting) with implementation of a procedure accounting for the difference between hybridization functions of the phys and aux systems. Such approach was implemented in [19] within the dual fermion technique. The technique employs an auxiliary fermionic degree of freedom—dual fermion—which provides a way to account for the difference in the form of diagrammatic expansion. This allows to get accurate results employing smaller number of auxiliary modes. However, additional procedure (superperturbation theory) is required on top of the QME solution (see [19] and references therein for technical details). We note in passing that the dual fermion approach relies on ability to evaluate multi-time correlation function in the aux system representing phys dynamics. Thus, proof of equivalence of multi-time correlation functions in phys and aux systems presented above is important for building theoretical foundation for the aux-DF method of [19].

Here we present a numerical simulation illustrating the equivalence of original unitary and Lindblad-type Markov treatment for the open quantum Fermi system. We note that the example is a simple illustration only and that realistic simulations will require more than two auxiliary modes and/or procedure to correct for approximate hybridization function to reproduce dynamics in the physical system. Such considerations have been done by us previously [19, 24].

We consider the Anderson impurity model (figure 3(a))

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{\sigma \in \{\uparrow ,\downarrow \}}{\varepsilon }_{0}{\hat{d}}_{\sigma }^{\dagger }{\hat{d}}_{\sigma }+U{\hat{n}}_{\uparrow }{\hat{n}}_{\downarrow }+\displaystyle \sum _{k\in L,R}\displaystyle \sum _{\sigma \in \{\uparrow ,\downarrow \}}\left({\varepsilon }_{k}{\hat{c}}_{k\sigma }^{\dagger }{\hat{c}}_{k\sigma }+{V}_{k}{\hat{d}}_{\sigma }^{\dagger }{\hat{c}}_{k\sigma }+{V}_{k}^{* }{\hat{c}}_{k\sigma }^{\dagger }{\hat{d}}_{\sigma }\right),\end{eqnarray} \tag{ 31 }$$

where ${\hat{n}}_{\sigma }={\hat{d}}_{\sigma }^{\dagger }{\hat{d}}_{\sigma }$. We calculate the system evolution after connecting initially empty site to baths at time t = 0. Parameters of the simulations are (numbers are in arbitrary units of energy E₀): ε₀ = 0 and U = 1. We assume

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{K}(E)={\gamma }_{K}\displaystyle \frac{{t}_{K}^{2}}{{\left(E-{\varepsilon }_{K}\right)}^{2}+{\left({\gamma }_{K}/2\right)}^{2}},\end{eqnarray} \tag{ 32 }$$

where ${{\rm{\Gamma }}}_{K}(E)\equiv 2\pi {\sum }_{k\in K}| {V}_{k}{| }^{2}\delta (E-{\varepsilon }_{k})$ is the electron escape rate into contact K (K = L, R), ε_L = ε_R = 0, γ_L = γ_R = 0.2, and t_L = t_R = 1.

Zoom In Zoom Out Reset image size

For simplicity, we consider infinite bias, so that auxiliary model with only two sites (figure 3(b)) is sufficient to reproduce dynamics in the physical system. After mapping, ε_L and ε_R become on-site energies of the auxiliary sites and γ_L and γ_R are taken as dissipation rates due to coupling to the L and R baths, respectively. In the auxiliary model we compare unitary evolution calculated within numerically exact td-DMRG [8, 9, 30, 31] with Lindblad QME results. Time is shown in units of ${t}_{0}={\hslash }/{E}_{0}$, currents in units of ${I}_{0}={E}_{0}/{\hslash }$, and ${\hslash }$ is assumed to be 1. Figure 4 shows level population, ${n}_{0}=\langle {\hat{n}}_{\sigma }\rangle$, as well as left, I_L, and right, I_R, currents in the system after quench. Close correspondence between the two numerical results is an illustration for exact analytical derivations presented in section 3.

Zoom In Zoom Out Reset image size

We consider an open quantum Fermi system S coupled to a number of external Fermi baths each at its own equilibrium (each bath has its own electrochemical potential μ_i and temperature T_i). The evolution of the model (system plus baths) is unitary. We show that reduced dynamics of the system S in the original unitary non-Markov model can be exactly reproduced by Markov non-unitary Lindblad-type evolution of an auxiliary system, which consists of the system S coupled to a number of auxiliary modes A which in turn are coupled to two Fermi baths L and R: one full (${\mu }_{L}\to +\infty$) and one empty (${\mu }_{R}\to -\infty$). The proof is performed in two steps: first we show that reduced S dynamics in the physical model is equivalent to reduced dynamics of S in the auxiliary model, when A degrees of freedom and the two baths are traced out; second, we show that reduced dynamics of S + A in the auxiliary model with unitary evolution of the model can be exactly reproduced by the Lindblad-type Markov evolution of S + A. The correspondence is shown to hold for reduced density matrix and for multi-time correlation functions defined on the Keldysh contour. Our study extends a recent work about Bose systems [26] to open Fermi systems and beyond only reduced density matrix consideration. Establishing the possibility of exact mapping of reduced unitary non-Markov dynamics to much simpler non-unitary Markov Lindbald-type treatment sets firm basis for auxiliary QME methods employed in, e.g. AMEA [17] or aux-DF [19] approaches. We note that in practical implementations improving the quality of mapping can be based on increasing number of A modes, as is done in advanced AMEA implementations [28], or by utilization of expansion in discrepancy between physical and auxiliary hybridization functions, as is done in the dual fermion formulation [29], or both. Scaling performance of the two approaches to mapping quality enhancement is a goal for future research.

We thank Max Sorantin for useful discussions. This material is based upon work supported by the National Science Foundation under grant CHE-1565939 (MG). EA gratefully acknowledges financial support by the Austrian Science Fund (FWF), grant No. P26508.

Here we prove that reduced density matrix of S + A in the aux model satisfies Markov Lindblad-type equation-of-motion (EOM), equation (20).

We start by considering unitary evolution of the aux model. Heisenberg EOM for bath annihilation operator ${\hat{c}}_{{Ck}}$ is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\hat{c}}_{{Ck}}(t)={\rm{i}}[{\hat{H}}^{{\rm{aux}}}(t),{\hat{c}}_{{Ck}}(t)]=-{\rm{i}}{\varepsilon }_{{Ck}}{\hat{c}}_{{Ck}}(t)-{\rm{i}}\displaystyle \sum _{m\in A}{\left({t}_{m}^{C}\right)}^{* }\,{\hat{a}}_{m}(t).\end{eqnarray} \tag{ A.1 }$$

Its formal solution is

$$\begin{eqnarray}&&{\hat{c}}_{{Ck}}(t)={{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Ck}}t}{\hat{c}}_{{Ck}}(0)-{\rm{i}}\displaystyle \sum _{m\in A}{\left({t}_{m}^{C}\right)}^{* }{\int }_{0}^{t}{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}{\varepsilon }_{{Ck}}(t-s)}{\hat{a}}_{m}(s).\end{eqnarray} \tag{ A.2 }$$

Thus, Heisenberg EOM for an arbitrary operator $\hat{O}$ on S + A can be written as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\hat{O}(t)={\rm{i}}[{\hat{H}}_{{SA}}(t),\hat{O}(t)]-{\rm{i}}\displaystyle \sum _{m\in A}\left\{\displaystyle \sum _{k\in R}\left[{t}_{m}^{R}{\left[\hat{O}(t),{\hat{a}}_{m}^{\dagger }\right]}_{\zeta }\left({{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Rk}}t}{\hat{c}}_{{Rk}}(0)-{\rm{i}}\displaystyle \sum _{m^{\prime} \in A}{\left({t}_{m^{\prime} }^{R}\right)}^{* }{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Rk}}(t-s)}{\hat{a}}_{m^{\prime} }(s)\right)\right.\right.\\ \quad \ +\ \left.\zeta {\left({t}_{m}^{R}\right)}^{* }\left({{\rm{e}}}^{{\rm{i}}{\varepsilon }_{{Rk}}t}{\hat{c}}_{{Rk}}^{\dagger }(0)+{\rm{i}}\displaystyle \sum _{m^{\prime} \in A}{t}_{m^{\prime} }^{R}{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}{\varepsilon }_{{Rk}}(t-s)}{\hat{a}}_{m^{\prime} }^{\dagger }(s)\right){\left[\hat{O}(t),{\hat{a}}_{m}(t)\right]}_{\zeta }\right]\\ \quad \ -\ \displaystyle \sum _{k\in L}\left[{\left({t}_{m}^{L}\right)}^{* }{\left[\hat{O}(t),{\hat{a}}_{m}(t)\right]}_{\zeta }\left({{\rm{e}}}^{{\rm{i}}{\varepsilon }_{{Lk}}t}{\hat{c}}_{{Lk}}^{\dagger }(0)+{\rm{i}}\displaystyle \sum _{m^{\prime} \in A}{t}_{m^{\prime} }^{L}{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}{\varepsilon }_{{Lk}}(t-s)}{\hat{a}}_{m^{\prime} }^{\dagger }(s)\right)\right.\\ \quad \ +\ \left.\left.\zeta {t}_{m}^{L}\left({{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Lk}}t}{\hat{c}}_{{Lk}}(0)-{\rm{i}}\displaystyle \sum _{m^{\prime} \in A}{\left({t}_{m^{\prime} }^{L}\right)}^{* }{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Lk}}(t-s)}{\hat{a}}_{m^{\prime} }(s)\right){\left[\hat{O}(t),{\hat{a}}_{m}^{\dagger }(t)\right]}_{\zeta }\right]\right\},\end{array}\end{eqnarray} \tag{ A.3 }$$

where ζ = ±1 if $\hat{O}$ contains even/odd number of fermion operators, and [,]_ζ is (anti)commutator for $\zeta =(-)1$.

For future reference we introduce

$$\begin{eqnarray}&&{\hat{c}}_{C}^{({\rm{in}})}(t)\equiv \displaystyle \frac{1}{\sqrt{2\pi {N}_{C}}}\displaystyle \sum _{k\in C}{{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Ck}}t}{\hat{c}}_{{Ck}}(0)\end{eqnarray} \tag{ A.4 }$$

which satisfies anti-commutation relations

$$\begin{eqnarray}&&\left\{{\hat{c}}_{{C}_{1}}^{({\rm{in}})}(t),{\hat{c}}_{{C}_{2}}^{({\rm{in}})\dagger }(s)\right\}={\delta }_{{C}_{1},{C}_{2}}\,\delta (t-s),\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&\left\{{\hat{c}}_{{C}_{1}}^{({\rm{in}})}(t),{\hat{c}}_{{C}_{2}}^{({\rm{in}})}(s)\right\}=\left\{{\hat{c}}_{{C}_{1}}^{({\rm{in}})\dagger }(t),{\hat{c}}_{{C}_{2}}^{({\rm{in}})\dagger }(s)\right\}=0.\end{eqnarray} \tag{ A.6 }$$

The delta in time is due to the fact that the contact density of states N_C is constant, and, thus

$$\begin{eqnarray}&&\displaystyle \sum _{k\in C}{{\rm{e}}}^{-{\rm{i}}{\varepsilon }_{{Ck}}t}\equiv \int {\rm{d}}\varepsilon \,{N}_{C}(\varepsilon ){{\rm{e}}}^{-{\rm{i}}\varepsilon t}=2\pi {N}_{C}\delta (t)\end{eqnarray} \tag{ A.7 }$$

is satisfied. Note also that

$$\begin{eqnarray}&&{\int }_{0}^{t}{\rm{d}}{s}\,\delta (t-s)f(s)=\displaystyle \frac{1}{2}f(t)\end{eqnarray} \tag{ A.8 }$$

holds for arbitrary function f(t).

Using (A.4), (A.7), and (A.8) in (A.3) leads to

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\hat{O}(t)={\rm{i}}[{\hat{H}}_{{SA}}(t),\hat{O}(t)]\\ \quad \ -\ {\rm{i}}\displaystyle \sum _{m\in A}\left\{{t}_{{Rj}}{\left[\hat{O}(t),{\hat{a}}_{m}^{\dagger }(t)\right]}_{\zeta }\,{\hat{c}}_{R}^{({\rm{in}})}(t)+\zeta {\left({t}_{m}^{R}\right)}^{* }{\hat{c}}_{R}^{({\rm{in}})\dagger }(t){\left[\hat{O}(t),{\hat{a}}_{m}(t)\right]}_{\zeta }\right.\\ \quad \ -\ \left.{\left({t}_{m}^{L}\right)}^{* }{\left[\hat{O}(t),{\hat{a}}_{m}(t)\right]}_{\zeta }{\hat{c}}_{L}^{({\rm{in}})\dagger }(t)-\zeta {t}_{m}^{L}{\hat{c}}_{L}^{({\rm{in}})}(t){\left[\hat{O}(t),{\hat{a}}_{m}^{\dagger }(t)\right]}_{\zeta }\right\}\\ \quad \ -\ \displaystyle \frac{1}{2}\displaystyle \sum _{{m}_{1},{m}_{2}\in A}\left\{{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}^{\dagger }(t)\right]}_{\zeta }\,{\hat{a}}_{{m}_{2}}(t)-\zeta {{\rm{\Gamma }}}_{{m}_{2}{m}_{1}}^{R}{\hat{a}}_{{m}_{2}}{\left(t\right)}^{\dagger }{\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}(t)\right]}_{\zeta }\right.\\ \quad \ -\ \left.\zeta {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\hat{a}}_{{m}_{2}}(t){\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}^{\dagger }(t)\right]}_{\zeta }+{{\rm{\Gamma }}}_{{m}_{2}{m}_{1}}^{L}{\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}(t)\right]}_{\zeta }\,{\hat{a}}_{{m}_{2}}^{\dagger }(t)\right\},\end{array}\end{eqnarray} \tag{ A.9 }$$

where we employed definition of the dissipation matrix, equation (22).

Next we are going to write EOM for expectation value of $\hat{O}$

$$\begin{eqnarray}&&\left\langle \hat{O}(t)\right\rangle \equiv {\rm{Tr}}\left[\hat{O}(t){\hat{\rho }}^{{\rm{aux}}}(0)\right]\end{eqnarray} \tag{ A.10 }$$

by averaging (A.9) with initial density operator of the aux model, equation (19). Because initially S + A is from the baths and because bath L is full and R is empty (see figure 1(b))

$$\begin{eqnarray}&&{\hat{c}}_{L}^{({\rm{in}})\dagger }(t){\hat{\rho }}_{L}={\hat{\rho }}_{L}\,{\hat{c}}_{L}^{({\rm{in}})}(t)={\hat{c}}_{R}^{({\rm{in}})}(t){\hat{\rho }}_{R}={\hat{\rho }}_{R}\,{\hat{c}}_{R}^{({\rm{in}})\dagger }(t)=0\end{eqnarray} \tag{ A.11 }$$

holds. Thus, second and third lines in (A.9) do not contribute, and EOM for the expectation value of $\hat{O}(t)$ is

$$\begin{eqnarray}&&\begin{array}{l}\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\hat{O}(t)\right\rangle ={\rm{Tr}}\left[{\hat{\rho }}^{{\rm{aux}}}(0){\rm{i}}[{\hat{H}}_{{SA}}(t),\hat{O}(t)\right]-\displaystyle \frac{1}{2}\\ \quad \ \times \ \displaystyle \sum _{{m}_{1},{m}_{2}\in A}{\rm{Tr}}\left[{\hat{\rho }}^{{\rm{aux}}}(0)\left\{{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}^{\dagger }(t)\right]}_{\zeta }\,{\hat{a}}_{{m}_{2}}(t)-\zeta {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\hat{a}}_{{m}_{1}}^{\dagger }(t){\left[\hat{O}(t),{\hat{a}}_{{m}_{2}}(t)\right]}_{\zeta }\right.\right.\\ \quad \ \left.\left.+\ {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\left[\hat{O}(t),{\hat{a}}_{{m}_{2}}(t)\right]}_{\zeta }\,{\hat{a}}_{{m}_{1}}^{\dagger }(t)-\zeta {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\hat{a}}_{{m}_{2}}(t){\left[\hat{O}(t),{\hat{a}}_{{m}_{1}}^{\dagger }(t)\right]}_{\zeta }\right\}\right].\end{array}\end{eqnarray} \tag{ A.12 }$$

Because $\hat{O}$ is arbitrary in S + A, after transforming to Schrödinger picture (A.12) can be rewritten as EOM for ${\hat{\rho }}^{{\rm{aux}}}(t)$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\hat{\rho }}^{{\rm{aux}}}(t) & = & -{\rm{i}}[{\hat{H}}_{{SA}}(t),{\hat{\rho }}^{{\rm{aux}}}(t)]\\ & & +\ \displaystyle \sum _{{m}_{1},{m}_{2}\in A}\left[{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}\left({\hat{a}}_{{m}_{2}}{\hat{\rho }}^{{\rm{aux}}}(t){\hat{a}}_{{m}_{1}}^{\dagger }-\displaystyle \frac{1}{2}\left\{{\hat{\rho }}^{{\rm{aux}}}(t),{\hat{a}}_{{m}_{1}}^{\dagger }{\hat{a}}_{{m}_{2}}\right\}\right)\right.\\ & & +\ \left.{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}\left({\hat{a}}_{{m}_{1}}^{\dagger }{\hat{\rho }}^{{\rm{aux}}}(t){\hat{a}}_{{m}_{2}}-\displaystyle \frac{1}{2}\left\{{\hat{\rho }}^{{\rm{aux}}}(t),{\hat{a}}_{{m}_{2}}{\hat{a}}_{{m}_{1}}^{\dagger }\right\}\right)\right].\end{array}\end{eqnarray} \tag{ A.13 }$$

Finally, because only operators in S + A subspace appear in the right side of (A.13), tracing out baths degrees of freedom leads to equation (20).

Here we prove that two-time correlation function of two arbitrary operators in S + A, $\langle {\hat{O}}_{1}({t}_{1}){\hat{O}}_{2}({t}_{2})\rangle$ (t₁ ≥ t₂ ≥ 0), equation (23), can be equivalently obtained from reduced Lindblad-type evolution in the S + A subspace of the aux model.

Let introduce t ≡ t₁−t₂ ≥ 0, then ${\hat{O}}_{1}({t}_{1}){\hat{O}}_{2}({t}_{2})={\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2})$ and using equation (A.9) we get

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2}) & = & \left\{{\rm{i}}[{\hat{H}}_{{SA}}(t+{t}_{2}),{\hat{O}}_{1}(t+{t}_{2})]\right.\\ & & -\ {\rm{i}}\displaystyle \sum _{m\in A}\left({t}_{m}^{R}{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{m}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{c}}_{R}^{({\rm{in}})}(t+{t}_{2})\right.\\ & & +\ {\zeta }_{1}{\left({t}_{m}^{R}\right)}^{* }{\hat{c}}_{R}^{({\rm{in}})\dagger }(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{m}(t+{t}_{2})\right]}_{{\zeta }_{1}}\\ & & -\ {\left({t}_{m}^{L}\right)}^{* }{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{m}(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{c}}_{L}^{({\rm{in}})\dagger }(t+{t}_{2})\\ & & -\ \left.{\zeta }_{1}{t}_{m}^{L}{\hat{c}}_{L}^{({\rm{in}})}(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{m}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}\right)\\ & & -\ \displaystyle \frac{1}{2}\displaystyle \sum _{{m}_{1},{m}_{2}\in A}\left({{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right.\\ & & -\ {\zeta }_{1}{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right]}_{{\zeta }_{1}}\\ & & +\ {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\\ & & \ \left.\left.-{\zeta }_{1}{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\hat{a}}_{{m}_{2}}(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}\right)\right\}{\hat{O}}_{2}({t}_{2}).\end{array}\end{eqnarray} \tag{ B.1 }$$

Note that for t > 0

$$\begin{eqnarray}&&{\left[{\hat{c}}_{C}^{({\rm{in}})\dagger }(t+{t}_{2}),{\hat{O}}_{2}({t}_{2})\right]}_{{\zeta }_{2}}={\left[{\hat{c}}_{C}^{({\rm{in}})}(t+{t}_{2}),{\hat{O}}_{2}({t}_{2})\right]}_{{\zeta }_{2}}=0.\end{eqnarray} \tag{ B.2 }$$

Indeed, because from equation (A.9) it is clear that ${\hat{O}}_{2}({t}_{2})$ depends on ${\hat{O}}_{2}(s)$ and ${\hat{c}}_{C}^{({\rm{in}})(\dagger )}(s)$ only at earlier times (s < t₂) and because equation (A.5) shows that ${\hat{c}}_{C}^{({\rm{in}})(\dagger )}(s)$ taken at different times anti-commute with each other, equation (B.2) holds.

Thus, while taking the expectation value of (B.1) with respect to ${\hat{\rho }}^{{\rm{aux}}}(0)$, equation (19), ${\hat{c}}_{L}^{({\rm{in}})\dagger }(t+{t}_{2})$ and ${\hat{c}}_{R}^{({\rm{in}})}(t+{t}_{2})$ can be moved over ${\hat{O}}_{2}({t}_{2})$ for any t > 0. So as in appendix A, terms with ${\hat{c}}_{C}^{({\rm{in}})(\dagger )}(t)$ in (B.1) again don't contribute (see equation (A.11)), and we get for t > 0

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left\langle {\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2})\right\rangle ={\rm{Tr}}\left[\left\{{\rm{i}}[{\hat{H}}_{{SA}}(t+{t}_{2}),{\hat{O}}_{1}(t+{t}_{2})]\right.\right.\\ \quad \ -\ \displaystyle \frac{1}{2}\displaystyle \sum _{{m}_{1},{m}_{2}\in A}\left({{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right.\\ \quad \ -\ {\zeta }_{1}{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{R}{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right]}_{{\zeta }_{1}}\\ \quad \ +\ {{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{2}}(t+{t}_{2})\right]}_{{\zeta }_{1}}{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\\ \quad \ \left.\left.\left.-\ {\zeta }_{1}{{\rm{\Gamma }}}_{{m}_{1}{m}_{2}}^{L}{\hat{a}}_{{m}_{2}}(t+{t}_{2}){\left[{\hat{O}}_{1}(t+{t}_{2}),{\hat{a}}_{{m}_{1}}^{\dagger }(t+{t}_{2})\right]}_{{\zeta }_{1}}\right)\right\}\times {\hat{O}}_{2}({t}_{2}){\hat{\rho }}^{{\rm{aux}}}(0)\right].\end{array}\end{eqnarray} \tag{ B.3 }$$

Rearranging evolution operators, equation (24), and separating traces over S + A and L + R yields

$$\begin{eqnarray}&&\left\langle {\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2})\right\rangle ={{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{1}\,{{\rm{Tr}}}_{{LR}}\left[{\hat{U}}^{{\rm{aux}}}(t+{t}_{2},0){\hat{O}}_{2}({t}_{2}){\hat{\rho }}^{{\rm{aux}}}(0){\hat{U}}^{{\rm{aux}}\dagger }(t+{t}_{2})\right]\right\},\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left\langle {\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2})\right\rangle ={{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{1}\,\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{\rm{Tr}}}_{{LR}}\left[{\hat{U}}^{{\rm{aux}}}(t+{t}_{2},0){\hat{O}}_{2}({t}_{2}){\hat{\rho }}^{{\rm{aux}}}(0){\hat{U}}^{{\rm{aux}}\dagger }(t+{t}_{2})\right]\right\}.\end{eqnarray} \tag{ B.5 }$$

So that (B.3) can be rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{1}\,\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{\rm{Tr}}}_{{LR}}[\ldots ]\right\} & = & {{\rm{Tr}}}_{{SA}}\left\{\left({{ \mathcal L }}_{{SA}}^{\dagger }(t){\hat{O}}_{1}\right){{\rm{Tr}}}_{{LR}}[\ldots ]\right\}\\ & \equiv & \ {{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{1}\,{{ \mathcal L }}_{{SA}}(t){{\rm{Tr}}}_{{LR}}[\ldots ]\right\},\end{array}\end{eqnarray} \tag{ B.6 }$$

where ${{ \mathcal L }}_{{SA}}^{\dagger }(t)$ is adjoint [32] of the Liouvillian ${{ \mathcal L }}_{{SA}}(t)$ defined in (20), and where ${{\rm{Tr}}}_{{LR}}[\ldots ]$ is used as a shorthand notation for the full expression in (B.4), (B.5).

Taking into account that ${\hat{O}}_{1}$ is an arbitrary operator, we get

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{\rm{Tr}}}_{{LR}}[\ldots ]={{ \mathcal L }}_{{SA}}\,{{\rm{Tr}}}_{{LR}}[\ldots ]\end{eqnarray} \tag{ B.7 }$$

which has solution

$$\begin{eqnarray}&&{{\rm{Tr}}}_{{LR}}[\ldots ](t)={{ \mathcal U }}_{{SA}}(t,0){{\rm{Tr}}}_{{LR}}[\ldots ](0)\equiv {{ \mathcal U }}_{{SA}}(t,0){\hat{O}}_{2}\,{\hat{\rho }}_{{SA}}({t}_{2}).\end{eqnarray} \tag{ B.8 }$$

Substituting (B.8) into (B.4) leads to

$$\begin{eqnarray}&&\left\langle {\hat{O}}_{1}(t+{t}_{2}){\hat{O}}_{2}({t}_{2})\right\rangle ={{\rm{Tr}}}_{{LR}}\left[{\hat{O}}_{1}\,{{ \mathcal U }}_{{SA}}({t}_{1},{t}_{2})\left({\hat{O}}_{2}\,{\hat{\rho }}_{{SA}}({t}_{2})\right)\right].\end{eqnarray} \tag{ B.9 }$$

This relation expresses two-time correlation function defined from unitary evolution of the aux model in terms of Lindblad-type evolution of S + A subspace of the aux model. Finally, we note that while we had restriction t > 0 in derivation of (B.3), the result is correct also for t = 0, as one can see by direct comparison of the two sides in (B.9). Equation (B.9) together with (20) leads to (25).

Similarly, for t₂ ≥ t₁ ≥ 0 one can prove that

$$\begin{eqnarray}&&\langle {\hat{O}}_{1}({t}_{1}){\hat{O}}_{2}({t}_{2})\rangle =\langle \langle I| {{ \mathcal O }}_{2}^{-}\,{{ \mathcal U }}_{{SA}}({t}_{1},{t}_{2}){{ \mathcal O }}_{1}^{+}\,{{ \mathcal U }}_{{SA}}({t}_{2},0)| {\rho }_{{SA}}(0)\rangle \rangle .\end{eqnarray} \tag{ B.10 }$$

Here we prove that multi-time correlation functions of arbitrary operators ${\hat{O}}_{i}$ in S + A of the aux model defined on the Keldysh contour

$$\begin{eqnarray}&&\left\langle {T}_{c}{\hat{O}}_{1}({\tau }_{1}){\hat{O}}_{2}({\tau }_{2})\ldots \,{\hat{O}}_{N}({\tau }_{N})\right\rangle ,\end{eqnarray} \tag{ C.1 }$$

can be evaluated from Markov Lindblad-type evolution in the S + A subspace. Here operators ${\hat{O}}_{i}$ are in the Heisenberg picture. Projections (one-the-contour time orderings) of multi-time correlation functions (C.1) will have the following form

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {\hat{B}}_{1}({s}_{1}){\hat{B}}_{2}({s}_{2})..{\hat{B}}_{m}({s}_{m}){\hat{C}}_{n}({t}_{n})...{\hat{C}}_{2}({t}_{2}){\hat{C}}_{1}({t}_{1})\right\rangle \\ \quad \ =\ {\rm{Tr}}\left[{\hat{C}}_{n}({t}_{n})...{\hat{C}}_{2}({t}_{2}){\hat{C}}_{1}({t}_{1}){\hat{\rho }}^{{\rm{aux}}}(0){\hat{B}}_{1}({s}_{1}){\hat{B}}_{2}({s}_{2})...{\hat{B}}_{m}({s}_{m})\right],\end{array}\end{eqnarray} \tag{ C.2 }$$

where ${\hat{B}}_{j}({s}_{j})$ and ${\hat{C}}_{i}({t}_{i})$ are used for operators ${\hat{O}}_{i}$ on the backward and forward branches of the contour, respectively (see figure 2) and where

$$\begin{eqnarray}&&\begin{array}{l}{t}_{n}\gt {t}_{n-1}\gt ...\gt {t}_{1}\geqslant 0\\ {s}_{m}\gt {s}_{m-1}\gt ...\gt {s}_{1}\geqslant 0.\end{array}\end{eqnarray} \tag{ C.3 }$$

Note, there is no ordering between the sets {t_i} and {s_j} (i ∈ {1, 2, ..., n} and j ∈ {1, 2, ..., m}).

Let denote the time-ordering of the set {t₁, t₂, ..., t_n, s₁, s₂, ..., s_m} by {θ₁, ..., θ_m + n}. So that

$$\begin{eqnarray}&&{\theta }_{m+n}\geqslant {\theta }_{m+n-1}\geqslant \ldots \geqslant {\theta }_{1}\geqslant 0.\end{eqnarray} \tag{ C.4 }$$

We want to prove that projections of multi-time correlation functions satisfy quantum regression theorem [32]

$$\begin{eqnarray}&&\begin{array}{l}\langle \langle I| {{ \mathcal O }}_{{\theta }_{m+n}}\,{{ \mathcal U }}_{{SA}}({\theta }_{m+n},{\theta }_{m+n-1}){{ \mathcal O }}_{{\theta }_{m+n-1}}\,{{ \mathcal U }}_{{SA}}({\theta }_{m+n-1},{\theta }_{m+n-2})\ldots \\ \quad \ \qquad \ \ldots {{ \mathcal O }}_{{\theta }_{1}}\,{{ \mathcal U }}_{{SA}}({\theta }_{1},0)| {\rho }_{{SA}}(0)\rangle \rangle ,\end{array}\end{eqnarray} \tag{ C.5 }$$

where ${{ \mathcal O }}_{{\theta }_{i}}$ is superoperator, equation (26), corresponding to operator $\hat{B}$ or $\hat{C}$ (backward or forward branch of the contour, respectively) at real time θ_i.

We prove (C.5) by mathematical induction. First, we note that equations (25) and (B.10) are special cases of equation (C.5) with m + n = 2. Suppose that for any combination (m, n) satisfying m + n = k, equation (C.5) holds. Now let consider (k + 1)-time correlation function

$$\begin{eqnarray}&&\left\langle {\hat{B}}_{1}({s}_{1}){\hat{B}}_{2}({s}_{2})\ldots {\hat{B}}_{m}({s}_{m}){\hat{O}}_{{\theta }_{k+1}}({\theta }_{k+1}){\hat{C}}_{n}({t}_{n})\ldots {\hat{C}}_{2}({t}_{2}){\hat{C}}_{1}({t}_{1})\right\rangle ,\end{eqnarray} \tag{ C.6 }$$

where θ_k + 1 > t_n > t_n−1 > ... > t₁ ≥ 0 and θ_k + 1 > s_m > s_m−1 > ... > s₁ ≥ 0. As previously, we time-order both sets

$$\begin{eqnarray}&&{\theta }_{k+1}\gt {\theta }_{k}\geqslant {\theta }_{k-1}\geqslant \ldots \geqslant {\theta }_{1}\geqslant 0,\end{eqnarray} \tag{ C.7 }$$

and take the derivative with respect to the latest time

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{\theta }_{k+1}}\left\langle {\hat{B}}_{1}({s}_{1}){\hat{B}}_{2}({s}_{2})\ldots {\hat{B}}_{m}({s}_{m}){\hat{O}}_{k+1}({\theta }_{k+1}){\hat{C}}_{n}({t}_{n})\ldots {\hat{C}}_{2}({t}_{2}){\hat{C}}_{1}({t}_{1})\right\rangle \\ \quad \ \equiv \ \displaystyle \frac{{\rm{d}}}{{\rm{d}}{\theta }_{k+1}}\langle \langle I| {{ \mathcal O }}_{{\theta }_{k+1}}\,{{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},{\theta }_{k}){{ \mathcal O }}_{{\theta }_{k}}\,{{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k},{\theta }_{k-1})\ldots {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{1},0)| {\rho }^{{\rm{aux}}}(0)\rangle \rangle \\ \quad \ =\ {{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{{\theta }_{k+1}}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{\theta }_{k+1}}\langle \langle {I}_{{LR}}| {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},{\theta }_{k}){{ \mathcal O }}_{{\theta }_{k}}\,{{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k},{\theta }_{k-1})\ldots {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{1},0)| {\rho }^{{\rm{aux}}}(0)\rangle {\rangle }_{{LR}}\right\}\\ \quad \ =\ {{\rm{Tr}}}_{{SA}}\left\{{\hat{O}}_{{\theta }_{k+1}}{{ \mathcal L }}_{{SA}}({\theta }_{k+1})\langle \langle {I}_{{LR}}| {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},{\theta }_{k}){{ \mathcal O }}_{{\theta }_{k}}\,{{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k},{\theta }_{k-1})\ldots {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{1},0)| {\rho }^{{\rm{aux}}}(0)\rangle {\rangle }_{{LR}}\right\},\end{array}\end{eqnarray} \tag{ C.8 }$$

where we followed the argument leading to (B.5) and (B.6) in appendix B. In (C.8) ${{ \mathcal U }}^{{\rm{aux}}}$ is the Liouville space analog of the Hilbert space evolution operator ${\hat{U}}^{{\rm{aux}}}$ defined in equation (24).

Taking into account that ${\hat{O}}_{{\theta }_{1}}$ is an arbitrary operator, we get

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{\theta }_{k+1}}\langle \langle {I}_{{LR}}| {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},0)\ldots | {\rho }^{{\rm{aux}}}(0)\rangle {\rangle }_{{LR}}\\ \quad \ =\ {{ \mathcal L }}_{{SA}}({\theta }_{k+1})\langle \langle {I}_{{LR}}| {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},0)\ldots | {\rho }^{{\rm{aux}}}(0)\rangle {\rangle }_{{LR}},\end{array}\end{eqnarray} \tag{ C.9 }$$

where $\langle \langle {I}_{{LR}}| {{ \mathcal U }}^{{\rm{aux}}}({\theta }_{k+1},0)\ldots | {\rho }^{{\rm{aux}}}(0)\rangle {\rangle }_{{LR}}$ is shorthand notation for the expression introduced in (C.8).

Solving (C.9) and utilizing quantum regression theorem for its initial condition, θ_k + 1 = θ_k, leads to

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {\hat{B}}_{1}({s}_{1}){\hat{B}}_{2}({s}_{2})\ldots {\hat{B}}_{m}({s}_{m}){\hat{O}}_{{\theta }_{k+1}}({\theta }_{k+1}){\hat{C}}_{n}({t}_{n})\ldots {\hat{C}}_{2}({t}_{2}){\hat{C}}_{1}({t}_{1})\right\rangle \\ \quad \ =\ \langle \langle I| {{ \mathcal O }}_{{\theta }_{k+1}}\,{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{k+1}},{t}_{{\theta }_{k}}){{ \mathcal O }}_{{\theta }_{k}}\,{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{k}},{t}_{{\theta }_{k-1}})\ldots {{ \mathcal O }}_{{\theta }_{1}}{{ \mathcal U }}_{{SA}}({t}_{{\theta }_{1}},0)| {\rho }_{{SA}}(0)\rangle \rangle \end{array}\end{eqnarray} \tag{ C.10 }$$

which is quantum regression theorem for (k + 1)-time correlation function. Thus, by induction we prove equation (30).