Markovian Treatment of non-Markovian Dynamics of Open Fermionic Systems

We show that an open fermionic system coupled to continuous environment with unitary system-environment evolution can be exactly mapped onto an auxiliary system consisting of the physical fermion system and a set of discrete fermionic modes subject to non-unitary Lindblad-type system-modes evolution in such a way that reduced dynamics of the fermionic system in the two cases are the same. Conditions for equivalence of reduced dynamics in the two systems are identified and a proof is presented. The study is extension of recent work on Bose systems [D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 120, 030402 (2018)] to open quantum Fermi systems and to multi-time correlation functions. Numerical simulations within generic junction model are presented for illustration.


Introduction
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][2][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][13][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][21][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][24][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.

Models
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 iŝ Hereˆ( ) H t S andĤ B (Bä{1, K, N}) are Hamiltonians of the system and baths.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 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.
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 (m  +¥ L ), the other (R) is completely empty (m  -¥ R ). The Hamiltonian of the total system iŝ whereĤ S is the same as in (1) andV SA their interaction with the system Hereˆ † a m (â m ) creates (annihilates) electron in the auxiliary mode m in A. Ĥ C represents continuum of states in contact Ĉˆˆ( andV AC couples auxiliary modes A to bath C (L or R) 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.

Non-Markov to Markov mapping
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  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 The parameters of the auxiliary model should then be chosen such that the total hybridization functions for the system are as close as possible to the corresponding hybridization functions,˜( ) S E r and˜( ) S E K , 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)˜( 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 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,r SA , which is defined by integrating out the baths degrees of freedom from the total density operatorˆ( ) r aux follows an unitary evolution with initial condition given by S+A decoupled from the bathŝ arbitrary. In appendix A we prove thatˆ( ) r t SA satisfies the following Markov Lindblad-type equation of motion whereˆ(  SA is the Liouvillian superoperator defined on the S+A subspace of the aux model and 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Ô 1 andÔ 2 in S+A we define two-time (t 1 t 2 0) correlation function aŝ HereÛ aux is the evolution operator in the aux system 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 subspacê ,0 0 .
Finally, we extend consideration to multi-time correlation functions of arbitrary operatorsÔ i (iä{1, K, N}) defined on the Keldysh contour (see figure 2) aŝ where τ i are the contour variables, T c is the contour ordering operator, and 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 Here P is number of Fermi interchanges in the permutation of operatorsÔ i by T c , θ i are indices of operatorsÔ i rearranged is such a way that > > >  ( q t i is real time corresponding to contour variable t q i ), q  i are the superoperators defined in (26), and  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

Numerical illustration
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))ˆˆˆˆ(ˆˆˆˆˆˆ) 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 ): ε 0 =0 and U=1. We assume is the electron escape rate into contact K (K=L, R), ε L =ε R =0, γ L =γ R =0.2, and t L =t R =1.
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 , and  is assumed to be 1. Figure 4 shows level population,= á ñ s n n 0 , 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.

Conclusions
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 (m  +¥ L ) and one empty (m  -¥ R ). 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.
Here we prove that reduced density matrix of S+A in the aux model satisfies Markov Lindblad-type equationof-motion (EOM), equation (20).
We start by considering unitary evolution of the aux model. Heisenberg EOM for bath annihilation operator c Ck is   For future reference we introducê The delta in time is due to the fact that the contact density of states N C is constant, and, thus  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).

Appendix B. Derivation of equation (25)
Here we prove that two-time correlation function of two arbitrary operators in S+A,   Note, there is no ordering between the sets {t i } and {s j } (iä{1, 2, K, n} and jä{1, 2, K, m}).

a t t O t t a t t O t t a t t a t t a t t O t t a t t O t
Let denote the time-ordering of the set {t 1 , t 2 , K, t n , s 1 , s 2 , K, s m } by {θ 1 , K, θ m+n }. where q  i is superoperator, equation (26), corresponding to operatorB orĈ (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)