# Markovian Treatment of non-Markovian Dynamics of Open Fermionic Systems

###### Abstract

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.

## I Introduction

Open nonequilibirum systems are at the forefront of experimental and theoretical research due to 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 Jiang et al. (2009); Khasminskaya et al. (2016); Gaita-Ariño et al. (2019). 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 Anders (2008); Schmitt and Anders (2010), flow equations Wegner (1994); Kehrein (2006), time-dependent density matrix renormalization group Schollwöck (2005, 2011), multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) Wang and Thoss (2009, 2018), and continuous time quantum Monte Carlo Cohen et al. (2015); Antipov et al. (2017); Ridley et al. (2018) 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) Anisimov and Izyumov (2010); Aoki et al. (2014). 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) Arrigoni et al. (2013); Dorda et al. (2014) introducing numerically inexpensive and pretty accurate solver for the nonequilibrium DMFT. Another example is recent formulation of the auxiliary dual-fermion method Chen et al. (2019). While the mappings appear to be very useful and accurate, only semi-quantitative arguments for possibility of the mapping were presented with main supporting evidence being benchmarking vs. numerically exact computational techniques. In particular, a justification for the mapping was argued in Refs. Schwarz et al. (2016); Dorda et al. (2017); Arrigoni and Dorda (2018) based upon the singular coupling derivation of the Lindblad equation. Still, the consideration is not rigorous.

Recently, a rigorous proof of non-Markov to Markov mapping for open Bose quantum systems was presented in the literature Tamascelli et al. (2018). It was shown that evolution of reduced density matrix in non-Markov system with unitary system-environment evolution can be equivalently presented by Markov evolution of an extended system (system plus modes of environment) under non-unitary (Lindblad-type) evolution. Here, we extend consideration of Ref. Tamascelli et al. (2018) 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 II we discuss non-Markov to Markov mapping procedure in Section III. Exact mathematical proofs are given in Appendices. Section IV presents numerical illustration of the mapping for a simple generic model of a junction. We conclude in Section V.

## Ii Models

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

(1) |

Here and () are Hamiltonians of the system and baths. introduces coupling of the system to bath . 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.

(2) | |||||

(3) |

where () and () create (annihilate) electron in level of the system and level of bath . In the model, dynamics of the system-plus-baths evolution is unitary. Below we call this model (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’ll consider is a model we shall call (auxiliary; see Fig. 1b). Here, the same system is coupled to a number of auxiliary modes , which in their turn are coupled to two baths. There are two Fermi baths in the configuration: one () is completely full (), the other () is completely empty (). The Hamiltonian of the system is

(4) |

where is the same as in (1), represents set of modes

(5) |

and their interaction with the system

(6) |

Here () creates (annihilates) electron in the auxiliary mode in .

represents continuum of states in contact

(7) |

with constant density of states

(8) |

and couples auxiliary modes to bath ( or )

(9) |

Dynamics of the whole configuration is unitary.

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

## Iii Non-Markov to Markov mapping

First, we are going to prove that with an appropriate choice of parameters of the dynamics of can be equivalently represented in the original model and auxiliary model , 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 Antipov et al. (2017). 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 and are matrices introducing the corresponding hybridization functions for bath of the physical problem (Fig. 1a).

(10) | |||||

(11) |

where are the Fourier transforms of retarded (Keldysh) projections of the free electron Green’s function . Then total hybridization functions for the system

(12) | |||||

(13) |

should be identical with the corresponding hybridization functions, and , of in the auxiliary model (Fig. 1b). The latter have contribution from full () and empty () baths, and from auxiliary modes ()

(14) | |||||

(15) |

where we assume modes initially in stationary state. Requirement of equivalence can be expressed as

(16) | |||||

(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. In principle, this fitting can be done in many different ways Dorda et al. (2017). For example, possibility of exact fitting of an arbitrary function with set of Lorentzians was discussed in Ref. Imamoglu (1994). 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, as is implemented in, e.g, AMEA Dorda et al. (2015), or by employing diagrammatic expansion related to the difference between true and fitted hybridization functions, as is realized in, e.g., dual fermion approach Jung et al. (2012), or both.

Now, when equivalence of reduced system () dynamics in and is established, we turn to consideration of evolution of the model. We will show that reduced dynamics derived from unitary evolution of the model can be exactly represented by non-unitary Lindblad-type evolution.

Following Ref. Tamascelli et al. (2018) we consider reduced density operator of in , , which is defined by integrating out baths degrees of freedom of the total density operator

(18) |

The latter follows unitary evolution with initial condition being decoupled from the baths

(19) |

where and .

In Appendix A we prove that satisfies Markov Lindblad-type equation of motion

(20) | |||

where

(21) |

is the Liouvillian superoperator defined on the subspace of the model and

(22) |

is the dissipation matrix due to coupling to contact .

Next we turn to multi-time correlation functions of operators in the subspace of the model. Following Ref. Tamascelli et al. (2018) we start consideration from two-time correlation function on real time axis. For arbitrary operators and in we define two-time () correlation function as

(23) | |||

Here is the evolution operator in the system

(24) |

and is the time-ordering operator. In B we show that (23) can be equivalently obtained from reduced Linblad-type evolution in the subspace

(25) | |||

Here is Liouville space bra representation of the Hilbert space identity operator, is the Liouville space superoperator corresponding to the Hilbert space operator (see Fig. 2)

(26) |

and is the Liouville space evolution superoperator

(27) |

Finally, we extend consideration to multi-time correlation functions of arbitrary operators () defined on the Keldysh contour (see Fig. 2) as

(28) | |||

where are the contour variables, is the contour ordering operator, and

(29) |

is the contour evolution operator. Note subscripts of operators in the right side of (28) indicate both type of the operator and its position on the contour. In C we prove that multi-time correlation functions (28) can be evaluated solely from Markov Lindblad-type evolution in subspace of the model

(30) | |||

Here is number of Fermi interchanges in the permutation of operators by , are indices of operators rearranged is such a way that ( is real time corresponding to contour variable ), are the superoperators defined in (26), and is the Liouville space evolution superoperator defined in (27).

## Iv Numerical illustration

Here we present numerical simulation illustrating equivalence of original unitary and Lindblad-type Markov treatment for the open quantum Fermi system. We consider Anderson model (Fig. 3a)

(31) | ||||

where . We calculate the system evolution after connecting initially empty site to baths at time . Parameters of the simulations are (numbers are in arbitrary units of energy ): and . We assume

(32) |

where is the electron escape rate into contact , , and .

For simplicity, we consider high bias, so that auxiliary model with only two sites (Fig. 3b) is sufficient to reproduce dynamics in the physical system. In the auxiliary model we compare unitary evolution calculated within numerically exact td-DMRG Schollwöck (2005, 2011); Bauer et al. (2011); Dolfi et al. (2014) with Lindblad QME results. Time is shown in units of , currents use , and is assumed to be . Figure 4 shows level population, , as well as left, , and right, , currents in the system after quench. Close correspondence between the two numerical results is an illustration for exact analytical derivations presented in Section III.

## V Conclusions

We consider open quantum Fermi system coupled to a number of external Fermi baths each at its own equilibrium (each bath has its own electrochemical potential and temperature ). Evolution of the model (system plus baths) is unitary. We show that reduced dynamics of the system 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 coupled to a number of auxiliary modes which in turn are coupled to two Fermi baths and : one full () and one empty (). The proof is performed in two steps: first we show that reduced dynamics in the physical model is equivalent to reduced dynamics of in the auxiliary model, when degrees of freedom and the two baths are traced out; second, we show that reduced dynamics of + in the auxiliary model with unitary evolution of the model can be exactly reproduced by the Lindblad-type Markov evolution of +. The correspondence is shown to hold for reduced density matrix and for multi-time correlation functions defined on the Keldysh contour. Our study is extension of recent work about Bose systems Tamascelli et al. (2018) to open Fermi systems and beyond only reduced density matrix consideration. Establishing possibility of exact mapping of reduced unitary non-Markov dynamics to much simpler non-unitary Markov Lindbald-type treatment sets firm basis for auxiliary quantum master equations (QME) methods employed in, e.g, AMEA Arrigoni et al. (2013) or aux-DF Chen et al. (2019) approaches. We note that in practical implementations improving quality of mapping can be based on increasing number of modes, as is done in advanced AMEA implementations Dorda et al. (2015), or by utilization of expansion in discrepancy between physical and auxiliary hybridization functions, as is done in the dual fermion formulation Jung et al. (2012), or both. Scaling performance of the two approaches to mapping quality enhancement is a goal for future research.

###### Acknowledgements.

We thank Max Sorantin for useful discussions. This material is based upon work supported by the National Science Foundation under grant CHE-1565939.## Appendix A Derivation of Eq. (Iii)

Here we prove that reduced density matrix of in the model satisfies Markov Lindblad-type equation-of-motion (EOM), Eq. (III).

We start by considering unitary evolution of the model. Heisenberg EOM for bath annihilation operator is

(33) | ||||

Its formal solution is

(34) |

Thus, Heisenberg EOM for an arbitrary operator on can be written as

(35) | |||

where if contains even/odd number of fermion operators, and is (anti)commutator for .

For future reference we introduce

(36) |

which satisfies usual commutation relations

(37) | |||

(38) |

Note that because contact density of states is constant

(39) |

is satisfied. Note also that

(40) |

holds for arbitrary function .

Using (36), (39), and (40) in (35) leads to

(41) | ||||

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

Next we are going to write EOM for expectation value of

(42) |

by averaging (A) with initial density operator of the model, Eq. (19). Because initially is from the baths and because bath is full and is empty (see Fig. 1b)

(43) |

holds. Thus, second and third lines in (A) do not contribute, and EOM for the expectation value of is

(44) | |||

Because is arbitrary in , after transforming to Schrödinger picture (A) can be rewritten as EOM for

(45) | ||||

Finally, because only operators in subspace appear in the right side of (A), tracing out baths degrees of freedom leads to Eq. (III).

## Appendix B Derivation of Eq. (25)

Here we prove that two-time correlation function of two arbitrary operators in , (), Eq. (23), can be equivalently obtained from reduced Lindblad-type evolution in the subspace of the model.