Nonequilibrium Dynamical Mean Field Theory: an auxiliary Quantum Master Equation approach

We introduce a versatile method to compute electronic steady state properties of strongly correlated extended quantum systems out of equilibrium. The approach is based on dynamical mean-field theory (DMFT), in which the original system is mapped onto an auxiliary non-equilibrium impurity problem imbedded in a Markovian environment. The steady state Green's function of the auxiliary system is solved by full diagonalization of the corresponding Lindblad equation. The approach can be regarded as the nontrivial extension of the exact-diagonalization based DMFT to the non-equilibrium case. As a first application, we consider an interacting Hubbard layer attached to two metallic leads and present results for the steady-state current and the non-equilibrium density of states.

In this respect, the theoretical description and understanding of strongly correlated quantum systems out of equilibrium poses an exciting challenge to modern theoretical physics. A widely used and successful method to treat strongly correlated lattice systems in equilibrium is dynamical mean-field-theory [16][17][18] (DMFT). The success of the method lies on the one hand in the nontrivial treatment of dynamical properties, and on the other hand in its applicability to a range of different problems, from solid states fermionic systems to ultracold bosonic atoms, as well as the possibility to combine it with realistic electronic structure methods [19]. Recently, DMFT has been extended to deal with time dependent nonequilibrium problems [20][21][22][23][24][25]. The extensions are based on the Kadanoff-Baym-Keldysh nonequilibrium Green's function approach [26][27][28][29].
DMFT relies on the solution of a correlated impurity problem, which constitutes the bottleneck of the approach. Several techniques have been adopted in the equilibrium case. Most of them have been applied, in a more or less approximate or limited way, to nonequilibrium DMFT as well, either in steady state or within full time dependence. Among them are iterated perturbation theory [20], numerical renormalization group [23], continuous time quantum monte carlo (CTQMC) [24,30], noncrossing approximation and beyond [31]. Additionally, exact DMFT solutions are available in certain limits [21,32,33].
In this Letter, we propose an approach that, in contrast to previous work, while directly accessing steadystate properties, features a solution of the DMFT impurity problem with controlled accuracy. This means that the accuracy can be directly estimated by comparing the exact and the effective bath hybridization functions (Fig. 2). Also, no often unreliable analytical continuation is required.
At the heart of the method lies a solution of the nonequilibrium DMFT impurity problem, which can be seen as a generalization of the exact diagonalization (ED) approach, widely used in the equilibrium case [16]. However, a crucial difference with respect to conventional DMFT-ED is the fact that here the effective impurity model describes an infinite system and, thus, displays a continuous spectrum.
In ED-based equilibrium DMFT [16] a certain number of noninteracting bath sites is introduced in order to fit the bath hybridization function required by the the selfconsistency condition. The maximum number of bath sites is limited by the exponential increase of the manybody Hilbert space. In equilibrium, the fit is carried out in imaginary (Matsubara) frequency space, where functions are smooth, in contrast to real frequency. There are several difficulties in this approach when trying to extend it to nonequilibrium steady states. (i) Due to the finite number of bath sites, a stationary solution of the impurity problem will always produce some equilibrium self energy. Besides the fact that this may be questionable, it is not clear which chemical potential or temperature should be used for the impurity problem. (ii) Secondly, due to the finite system, the bath spectrum is discontinuous, so that a fit in real frequencies becomes problematic. Unfortunately, there is no Matsubara Green's function in nonequilibrium, so that this poses a serious problem.
The alternative presented in this Letter consists in replacing the DMFT impurity Hamiltonian with an effective one which is solvable by ED but at the same time describes a truly infinite system. This is obtained by connecting the interacting impurity to a moderate number of bath sites which, in turn, are attached to Markovian reservoirs, see below for details. The exact bath spectral function is smoothly obtained in the (ideal) limit of an infinite number of bath sites.
The action of such Markovian baths on the reduced density matrix of the system consisting of the other bath sites and of the impurity, is described by the Lindblad quantum Master equation [50]. The latter can be readily solved by diagonalizing the Lindbladian within the many-body "super-Fock" space of reduced density matrices of the system. Its solution determines both the retarded and Keldysh impurity Green's function, as well as the self-energy. The latter is used in the DMFT loop to obtain the new bath hybridization function, which is fitted by new bath parameters.
In order to illustrate the approach, we apply it to a simple model describing a heterojunction consisting of a correlated interface (c) sandwiched between two metallic leads α = l, r (see Fig. 1). Experimentally, such a setup has been recently explored where the correlated layer was realized by a V 2 O 3 microfilm that is coupled to Au leads [51].
Before a certain time τ < τ 0 the three regions c, l, r are disconnected and in equilibrium at different chemical potentials µ c , µ l , and µ r , respectively. This amounts to applying a bias voltage Φ = µ l − µ r between the leads. The central region, lying on the x − y plane, is described by a single-band Hubbard layer on a square lattice with onsite interaction U , onsite energy ε c , and nearest-neighbor hopping t. The leads consist of two half-infinite cubic lattices described by a nearest-neighbor noninteracting tight-binding model with hopping which we take as unit of the energy t = 1, and on-site energies ε α . We restrict for simplicity to the particle-hole symmetric case for which ε c = −U/2, µ r = −µ l , and ε r = −ε l . Finally, we take µ α = ε α , which corresponds to having the same electron densities in the two leads. A related nonequilibrium model has been treated in DMFT within perturbative impurity solvers in Refs. [25,31].
Starting at τ = τ 0 , a nearest-neighbor hopping v is switched on between the central region and the leads. After a sufficiently long time, a steady-state is reached, provided no trapped surface states occur. Nonequilibrium properties, in general, and nonlinear transport in particular can quite generally be addressed in the framework of the Keldysh Green's function approach [26][27][28][29]52]. Here, we adopt the notation (see e.g. Ref. [29]) where the (underlined) Keldysh Green's function is a 2 × 2 matrix containing the retarded (G R ), advanced (G A ), and Keldysh (G K ) components. In steady state, these depend on a single frequency only. The system is translation invariant in the direction parallel to the layer, which we denote as . Accordingly we can write Dyson's equation for the layer (c) Green's function G(k , ω), k being the momentum [52] as Here, Σ(k , ω) is the self-energy, g 0 is the v = 0, U = 0 layer Green's function, and g α (k , ω) are the v = 0 leads Green's functions on the first lead layers. Their retarded and Keldysh components are readily obtained analytically in terms of the Green's function of a halfinfinite tight-binding chain. Within DMFT, one approximates the self energy by a local, i.e. k -independent Σ(ω), which is determined by solving a quantum impurity model with the same interaction U embedded in a self-consistently determined bath [16]. The latter is completely specified by the bath hybridization function ∆(ω), which is determined selfconsistently by requiring that the Green's function of the impurity G IMP (ω) = G 0 (ω) −1 − Σ(ω) −1 be equal to the local Green's function of the layer (cf. [16,21,25]) where G(k , ω) is given by (1) with Σ(k , ω) = Σ(ω) as obtained by the solution of the impurity problem. Here, G 0 is the "Weiss" bare Green's function of the impurity model defined as The solution of the impurity problem is in fact the DMFT bottleneck. The usual (equilibrium) DMFT ED procedure consists in approximating the effect of the total bath hybridization function ∆(ω) by an "effective" bath with a finite number N b of bath sites. Quite generally, in equilibrium one carries out some fit to the bath hybridization function in Matsubara space. As discussed above, this is not appropriate in nonequilibrium steady state. In this Letter, we present a different approach: In additional to a certain (even) number N b of bath sites, which we more conveniently connect to the impurity in the form of two chain segments, we include two Markovian baths, which represent a par-ticle reservoir and sink, respectively. Their role is to compensate for the "missing" part of the infinite chain which would be necessary to exactly reproduce the desired ∆(ω). The bath parameters, i.e. hopping and on-site energies of the bath sites, as well as the Lindblad coefficients (see below) of the Markovian baths, are then fitted to ∆(ω). More specifically, we minimize the cost function dω while ∆ eff is the bath hybridization function produced by the effective bath (bath sites+Markovian baths). An important aspect is that, although the (outermost) baths are Markovian, their effect on the impurity site is non Markovian due to the presence of the intermediate bath sites. This can be seen, for example, in the spectrum of ∆ eff in Fig. 2, which in the Markovian case would be a constant. Furthermore, upon increasing the number N b of intermediate bath sites, the effect of the Markovian bath becomes weaker and one is expected to approach the exact result ∆ eff (ω) = ∆(ω) for large N b . We now specify the effective bath more in detail. This consists of an Hamiltonian for the "system" (a chain of impurity+bath sites) in usual notation, where 0 is the impurity site, and n = −1, · · · , −l (n = 1, · · · , l) are left (right) bath sites.
Here, E 0,0 = ε c , and the bath energies with real, symmetric, and positive definite Lindblad matrices Γ (i) n,m . At first sight, one would assume a "source" bath attached to the leftmost and a "sink" bath to the rightmost site. However, the accuracy improves considerably if one allows all Γ (i) n,m to be nonzero parameters to fit ∆.
As discussed in detail in Refs. [54][55][56][57], the opensystem problem describing the effective bath can be mapped onto a superhamiltonian acting on a superfermion space with twice as many degrees of freedom (i.e.,  "orbitals"). This many-body superhamiltonian, corresponding to iL, which is non-hermitian, can be diagonalized by conventional methods within the super-Hilbert space. Quite generally, L has a unique eigenvector with eigenvalue 0, which corresponds to the steady-state density matrix ρ SS . All other eigenvalues have a negative real part, corresponding to decaying terms. With the same formalism, and exploiting the quantum regression theorem [58], one can evaluate correlation functions C AB (τ ) = trA(τ )B(0)ρ SS of any pair of system operators A and B, and thus the required impurity self energy Σ(ω). [59] The noninteracting Green's function for the effective system+bath is necessary in order to extract Σ(ω) and the bath hybridization function ∆ eff . This can be easily obtained [60] by observing that the Markovian baths can be exactly represented by two noninteracting fermionic baths in the wide-band limit with chemical potentials ±∞. By taking into account the relation between the matrices Γ and the parameters of this bath [61], one obtains for the noninteracting system Green's function (boldface object represent matrices with indices corresponding to system sites n): The DMFT self-consistency loop consists in (i) starting from some initial values of the variational parameters E n,m and Γ (i) n,m , (ii) solving the impurity problem via the approach described above and determining Σ, (iii) evaluating G LOC (ω) and ∆(ω) (from (2)) (iv) determining new values of the parameters E n,m and Γ (i) n,m by minimizing the cost function, and finally (v) using these new parameters to repeat the procedure from (i) until the parameter values converge. Of course there is an intrinsic inaccuracy which, for a fixed number of bath sites cannot be reduced, and is due to the error in the fit of ∆(ω) by a finite number of parameters. In principle, this can be systematically improved by increasing the number of bath sites. Of course, this is limited by the exponential increase of the Hilbert space, which, in this case, is even faster due to the fact that the number of degrees of freedom of the superfermion space is twice the one of the fermion space, so this makes the effort more difficult than in ordinary ED for the same number of bath sites. One should, however, observe that the number of fit parameters is larger than in the case with "simple" ED without Markovian bath [62].
We have used here an effective bath containing N b = 2 sites. Still, by taking all possible parameters into account, i.e., allowing, for example, all matrix elements of the Γ matrices to be nonzero (within the constraints imposed by symmetries), this gives 8 independent fitting parameters. In "simple" ED, one would have only 2 in the particle-hole symmetric case [63]. We take model parameters [64] v 2 = 0.1, t = 1, and the leads are fixed to zero temperature. In Fig. 2, we show the result of the fit to the bath hybridization function for two values of U and bias voltage Φ = 2. As one can see, already with this small number of bath sites, the fit is quite good. Moreover, the structure of ∆ eff is clearly non-Markovian, as expected. For comparison, we also plot the result of the fit to the bath hybridization function obtained with N b = 4 bath sites. This shows a considerable improvement. The quality of the fit can be also inferred by directly plotting the local and the impurity Green's function in Fig. 3 for two different values of U and Φ = 2.
The bias voltage induces a steady-state current. The expression for the current density j (current per square plaquette) is obtained straightforwardly in terms of the layer Green's function and the v = 0 lead Green's func- tions (see, e.g. [52]) as: Re G R (k , ω)g K l (k , ω) This is plotted in Fig. 4 as a function of bias voltage. The current, as expected, decreases with increasing U for smaller biases. At larger Φ the behavior is opposite, since j extends over a range of voltages which increases with increasing U . While at U = 0, a particle going through the interface conserves k , and thus the current goes to zero at a bias voltage equal to the one-dimensional (zdirection) bandwidth, the scattering at nonzero U mixes k and thus broadens the bandwidth of possible final states. In Fig. 4 we plot the scaled current j/v 2 as a function of bias voltage for different values of U , as the conductance is expected to behave as v 2 in a conductor, while it is suppressed (∝ v 4 ) in a gapped system. The crossing of the curves around U ∼ 4 is a signal of the nearby equilibrium metal-insulator transition. However, this should be seen only as indicative, as the curves are taken at a relatively high bias.
In conclusion, we have introduced a versatile method to deal with strongly correlated systems out of equilibrium within dynamical mean-field theory. The DMFT self-consistent bath is approximated by an effective one consisting of a small number of sites coupled to a Markovian bath environment. The steady state and Green's function of the effective system is solved by ED of the corresponding Lindblad equation. The approach is particularly appropriate to deal directly with the steady state, without the need to consider full time evolution. Nevertheless, it should be straightforward, although computationally more demanding, to deal with time-dependent problems, e.g., to describe pump-probe processes.
The accuracy of the effective bath to reproduce the DMFT one obviously depends on the number of bath sites, which is limited by the exponential increase of the "super"-Hilbert space. Improvements can possibly go along solving the Lindblad problem in the "smaller" ordinary fermion space in combination with quantum trajectory methods [65][66][67], and/or by density matrix renormalization group approaches [68][69][70].
The approach illustrated here for a simple but experimentally relevant [51] model can be extended straightforwardly to a number of other physically relevant systems, including multi layer semiconducting heterostructures, ultracold atoms and correlated coupled-cavity arrays featuring driving and dissipation, molecular contacts, and can be used to study nonequilibrium quantum phase transitions in these systems. Extension to a nonlocal self energy, as in cluster DMFT, or in nonequilibrium variational/perturbative cluster approaches [71][72][73] is also an interesting development.
Conservation laws, as usual, reduce the size of the many-body Hilbert space one has to diagonalize. For additive conserved quantities, for example n σ (σ =↑, ↓), the density matrix "state" is characterized by the constraint n σ − n σ = 0. In general, one can define "sectors" with a given n σ − n σ , which can be separately diagonalized. For the Green's functions, these are characterized by n σ − n σ = ±1

C. Green's function
To evaluate the Green's function, one needs two time evolutions, i.e. expectation values of the form [58] G(B(t), A(t )) ≡ tr U B(t)A(t )ρ U = trB A t (t − t ) (6) for t > t , where ρ U is the density matrix of the universe (in contrast to ρ, which is the reduced density matrix of the "system"), and tr U = tr ⊗ tr E is the trace over the universe degrees of freedom, whereby tr is the one over the system, and tr E the one over the environment. Here, and H U is the total Hamiltonian of the universe.
The Quantum regression theorem [58] states that under the same assumptions for which (4) holds, one has (for s > 0) d ds A t (s) = LA t (s) .
In terms of eigenvalues and eigenvectors, this gives αL|A t (s) = e sLα αL|A t (0) .
For example, the "greater" Green's function iG >n,m (t, t ) ≡ tr U c n (t)c † m (t )ρ U = I | c n c † mt (t−t ) | I is obtained from (6) with B = c n , and A = c † m . Inserting the identity α | αR αL | and taking the steady state | ρ(t = +∞) = | α 0 R this gives iG >n,m (t + s, t ) = = α e (t−t )Lα α 0 L | c n | αR αL | c † m | α 0 R for t > t . In a similar way, one obtains the result for t < t and for the "lesser" Green's function. These are combined (and Fourier transformed), to obtain the retarded (G R (ω)), advanced G A (ω), and Keldysh G K (ω) components of the Green's function of the effective bath problem evaluated at the impurity site (n = m = 0).