Dual-Fermion approach to Non-equilibrium strongly correlated problems.
We present a generalization of the recently developed dual fermion approach introduced for correlated lattices to non-equilibrium problems. In its local limit, the approach has been used to devise an efficient impurity solver, the superperturbation solver for the Anderson impurity model (AIM). Here we show that the general dual perturbation theory can be formulated on the Keldysh contour. Starting from a reference Hamiltonian system, in which the time-dependent solution is found by exact diagonalization, we make a dual perturbation expansion in order to account for the relaxation effects from the fermionic bath. Simple test results for closed as well as open quantum systems in a fermionic bath are presented.
In the last years there has been an evergrowing experimental as well as theoretical interest to non-equilibrium physics in strongly correlated systems. Recent experiments Wall et al. (2009); Tobey et al. (2008) in the field of time resolved measurements have set tasks for theory that can not be fulfilled by the existing methods of analytical or numerical approaches to the considered systems. Apart from the time-dependent DMRGWhite and Feiguin (2004) and NRGAnders and Schiller (2005) methods, where a physical system is approximated by a finite one-dimensional chain, all other theoretical approaches are perturbative in a general sense, i.e. one subdivides the Hamiltonian of the system in two parts , solves the part exactly and takes the rest into account perturbatively.
Historically the first calculations of the transport through a quantum dot were the linear response investigation of the conductance by LandauerLandauer (1970). The further development of this kind of theorySchoeller and Schön (1994) lead to a series of strong coupling works. In those the perturbation expansion is done in powers of the tunneling amplitude (up to 2nd or 4th order) and the correlated part of the system is considered as an multi-orbital impurity. More sophisticated studies are based on the Keldysh technique. In the work of Wingreen and Meir Wingreen and Meir (1994) a general formula for the time evolution of a system as a function of the dot Green’s function and the tunneling is obtained and then the limit of weak symmetric tunneling is considered.
Another way to proceed is the weak-coupling approach, which starts from the transport through a non-interacting dot and applies Keldysh diagrammatics in terms of the Coulomb interaction. Analytical studies have been performed up to second order of the perturbation series and GW approximationThygesen (2008). Recent numerical studies based on a generalization of the continuous time Quantum Monte Carlo scheme (CT-QMC)Rubtsov et al. (2005) to the Keldysh contour Werner et al. (2010) are, in principle, also a natural generalization of the weak-coupling perturbation theory, where diagrams are sampled stochastically. The CT-QMC approaches are, though in principle exact, limited to rather short time scales, since the calculations are strongly affected by a dynamical phase problem in the non-equilibrium caseWerner et al. (2010).
Perturbative approaches all suffer from a common general drawback: they work well in some particular limit of parameter values, in which the system described by is a good starting point, and fail in the opposite limit. This is due to the fact that the unperturbed system in general is too simple to provide the basis of the perturbation theory in opposite limits.
One may instead consider the perturbation of a more complex system, simple enough to be solvable exactly and still sophisticated enough to include all important features of the considered full system so that all sources of complexity are possibly treated on equal footing. This task is in general non-trivial, as the systems that one would normally wish to take as the non-perturbed ones are strongly correlated and thus conventional diagrammatic techniques do not work.
The recently developed Rubtsov et al. (2008) dual-fermion technique was originally intended to solve lattice problems beyond the well-established DMFT approximation, i.e. to incorporate a -dependent into self-energies. One can consider it as an optimal perturbation around the DMFT starting point. In this technique, problems are solved in the above sense, i.e. a simple albeit nontrivial (reference) system is solved exactly, and the dual perturbation theory accounts for the difference between the original and the reference system. Therefore this technique has been called a superperturbation. By ”exactly” we mean finding the Green’s function and, in principle, all higher cumulants. In later works Hafermann et al. (2009a, b) it has been shown that the dual perturbation theory is indeed convergent in opposite limits and provides a sensible interpolation between these limits. For the case of the AIM, the dual perturbation theory has been to shown to pass into standard perturbation theory for the case of strong hybridization and weak interaction , while it reproduces the hybridization expansion in the opposite limit of weak hybridization and strong interaction. We expect a corresponding behavior for the approach considered here.
In this work we generalize the dual fermion approach to models out of equilibrium and present an implementation for the experimentally important case of a hopping quench. We use the path integral approach on the Keldysh contour to formulate the theory and to introduce a proper discretization scheme. We also present a simplified version of this method that is much less time- and computational power demanding but still proves to give qualitatively and sometimes even quantitatively correct results. The basic idea of the method is to generalize the previously discussed superperturbation impurity solver for the AIM. One replaces the continuous conduction electron bath by a small number of discrete bath levels (the reference system), solves this finite system by exact diagonalization (ED) and then accounts for the difference of hybridizations through the dual perturbation. In the non-equilibrium technique, the idea remains the same, only that now we work on the Keldysh contour and all the relevant correlation functions depend on the time variables themselves and not on time differences as, of course, in the non-stationary situation there is no time translation invariance.
The remainder of this work is organized as follows: In section II we derive the discretized action of a general system on the Keldysh contour and perform the dual decomposition. In section III we give the details of the implementation of the algorithm and present some benchmark results, comparing them to other methods and discussing them. In the last section we summarize our work.
Ii Description of the method
ii.1 Path integral formalism
It is well-known Keldysh (1964, 1965) that an accurate closed diagrammatic description of a non-equilibrium process is possible only on the so-called Keldysh contour Fig. 1. Here one propagates the state of the system first in the direction of increasing time along the time axis and then exactly backwards to the initial state. The reason for this is that in order to calculate averages of operators one has to propagate back to the initial state since the final state of the system, unlike in the zero-temperature ground state technique, is a priori unknown. Our goal will thus be to introduce the dual transformation for the Keldysh action.
In equilibrium quantum mechanics every observable is associated with a Hermitian operator . Its expectation value is given by , where is the density matrix of the system governed by the Hamiltonian . As long commutes with the Hamiltonian the expectation value of will not have any time dependence.
In the following we consider the situation in which the system is in equilibrium for times smaller than and is then perturbed by a sudden switch of some not specified internal parameter at . The theory at hand is in principle able to treat general time dependent perturbations, but this requires a significant numerical effort. For this case the expectation value of is given by the average of in the Heisenberg picture traced over the initial density matrix :
is the evolution operator of the system. If is the full, time dependent Hamiltonian of the system, the evolution operator is given by:
In the latter expression is the time ordering operator which reshuffles the operators into chronological order, with earlier times to the right. The anti-chronological time ordering operator rearranges later times to the right. Reading the time arguments of Eq. (1) from left to right, one sees that the time evolution of the observable starts at , propagates to t, and goes back again to and hence can be depicted by the Keldysh contour shown in Fig. 1. Introducing a time-ordering operator which arranges operators according to their position on the contour, the time-evolution operator can be written in the following form:
where the integral is taken along the contour. Considering the last part of Eq. (1), the time evolution of the system can be viewed as an initial value problem. At time the system is prepared by an initial density matrix and then the time evolution of the system is governed by Eq. (3). In cases where the system is correlated before time and the initial density matrix is unknown (i.e. cannot be written in the form , with some one-particle operator ), it is necessary to extend the contour in Fig. 1 to imaginary times as described in referenceWagner (1991). This can be done straightforwardly, but for simplicity we restrict ourselves to the case of non-correlated .
We start by deriving the Keldysh action following KamenevKamenev and Levchenko (2009) as it is most appropriate for our purposes. First we omit the interaction part and consider free electrons governed by the Hamiltonian , where small Greek letters enumerate the single-particle states.
The Keldysh partition function can formally be defined as
which is a complex way to express through the evolution operator along the Keldysh contour. The next step is to discretize the contour with points and insert the identity of the overcomplete fermion coherent-state basis Negele and Orland (1998)
at each discretization point . Here stands for coherent states corresponding to eigenvalues , small Latin letters enumerate the discretization points. We obtain
with given by:
In the latter expression the upper sign applies for and otherwise.
which implies that is non-correlated. in (II.1) is thus given by if the initial density matrix can be expressed as .
The expression for consists of three major parts: the first one describes the discrete time evolution backward in time along the lower branch of the contour, the second one the evolution forward and the third one represents the boundary condition of the fermionic states.
For reasons that will become apparent below, we will not take the limit , but seek for a suitable mapping of the Keldysh formalism on a time grid. Thus the above matrix form of the action is appropriate for us. The indices of the latter correspond to discretization points on the contour. It is not necessary to explicitly distinguish between the upper and lower branch of the contour and introduce the conventional Keldysh indices. The inverse Green’s function for two coupled non-interacting sites is shown in example 1.
Including of the interaction is straightforward, it results in terms of the type in the action (which is valid in the limit ). To finish of the preparation we notice that as baths are taken to be uncorrelated one can perform the Gaussian integration over the bath fields to obtain the action depending only on the impurity fields. This leads to the well-known hybridization function which in our case depends on two discrete time indices. Explicit calculation results in:
Here and must be considered as supermatrices in time indices and bath degrees of freedom. The structure of a single block of can be seen from example 1. While in most works the grid is chosen in such a way that all points are equidistant and the number of points on the upper and lower contour is equal (see, e.g., Ref. Freericks (2008)), we here need to choose the dimension of the time block to be odd due to the structure of the dual perturbation theory (see below).
ii.2 Dual perturbation theory on the Keldysh contour
In this section we generalize the superperturbation approach for the AIMHafermann et al. (2009b) to the case of non-equilibrium systems. We start from the path integral formulation of the partition function:
with the following discrete-time expression for the impurity action :
Here is the Green function of the isolated impurity and is the hybridization function, which describes the coupling of the impurity to a continuum of bath states.
To stress that we are working on a time grid, we write the indices of the time-dependent matrices explicitly. Indices for orbital and spin degrees of freedom have been summarized in Latin letters. is some local interaction.
We introduce a reference system with the same interaction part as the original system, albeit with coupling to a set of discrete bath levels described by a hybridization function . For a small number of bath states, this system can be solved exactly using ED.
The action of the original system can be expressed in terms of the reference system and a correction. To this end, the hybridization of the reference system is added and subtracted:
One has to be aware that in the non-equilibrium case the action representation always involves boundary conditions that enter the discrete matrices as shown in Example 1. When adding and subtracting we assume that both systems have the same contribution on the impurity. The other parts of the density matrix are less important, because the difference in these terms will be treated perturbatively by doing an expansion in . In the equilibrium case it was not necessary to discuss the boundary conditions of the action, because they were automatically fulfilled when working with Matsubara frequencies.
Now dual fermions (represented by and ) are introduced through a Gaussian identity for Grassmann variables. The transformation can be written in the following form:
with the following definitions for and :
being the Green function of the reference problem. After some straightforward algebra the partition function can be brought into the following form
The action can be split into three parts: Two parts, which contain either -fermions or dual variables, and a part that describes the coupling between them:
In the last equation we have combined temporal, orbital and spin indices into numbers. To integrate out the -fermion part, the following defining equation for the dual interaction potential is introduced:
Expanding the left hand side in and integrated over the -fermion fields produces the correlation functions of the reference system. The dual potential is evaluated by expanding the right hand side and comparing the expressions on both sides by order. Up to fourth order in -fermion fields we obtain:
where is the complete two-particle vertex of the reference problem:
being the full two-particle Green functions of the reference problem and
its disconnected part. The dual action can now be written as follows:
with the following definition of the bare dual Green’s function:
Now a usual diagrammatic expansion in powers of can be performed. The lowest order contribution to the dual self-energy is thus
The lowest order approximation to the Keldysh Green function is depicted diagrammatically in Fig. 2. After calculating an approximation to the dual propagator the result can be transformed back to c-fermions using the following exact relation:
Note that the factors in the definition of at each vertex of a diagram cancel with the corresponding factors in the expression for .
The above formulae require the inversion of the impurity Green function and hybridization. Therefore we employ a time discretization (note that, of course, the matrices are not diagonalized by Fourier transform).
Furthermore, there are some important side remarks for the numerical calculations. The first one involves the structure of the time grid. If it is chosen such that equally many equidistant points lie on the upper and lower branches of the contour and the total number of points is even, the time step between the last point on the upper contour and the next one on the lower contour is zero. This causes a vanishing matrix element in and leads to a singular matrix for . This can lead to a break down of the dual theory, because the hybridization has to be inverted. This problem can be cured if an additional point at the tip of the contour is introduced, which neither belongs to the upper nor to the lower contour. This additional point has no other consequences and can be treated without any further problems. That is the reason why unlike in the work Kamenev and Levchenko (2009), the dimension of the time block in Example 1 is odd. Another issue in the inversion procedure can occur, if the off-diagonal elements and of the density matrix vanish. This happens, if the impurity is totally decoupled from the bath for times smaller than . In that case in Eq. (10) becomes also ill conditioned. This problem can be cured by introducing an infinitesimally small hopping parameter for times smaller than .
ii.3 Calculation of one- and two-particle Green’s functions in exact diagonalization
In the last section it has become clear that the key ingredients to construct the superperturbation theory are the one and two particle Green’s function of the reference system. In the following we will give the formulae of those quantities in the Lehman representation. We start with the single-particle Green’s function. This quantity depends on two times and two indices:
To derive the spectral representation, identity matrices are inserted between the operators, and the time-evolution operator is written in its diagonal form:
and are times on the Keldysh contour, which means that the -functions are defined as follows:
where the relation symbols order the times along the contour.
In the above, denote the exact eigenstates of the reference problem, whereas the states account for the initial state of the system described by the density matrix . As mentioned above, we consider only such cases where the initial density matrix can be given as , where is some effective temperature and is a one-particle operator. Thus, states are the eigenstates of this auxiliary operator and are the corresponding eigenvalues.
The formal definition of the two particle Green’s function is given by the following expression:
Here it was necessary to introduce abbreviations for the operators, in order to rewrite the time-ordered product as a sum over all possible permutations multiplied by a -function depending on the four time arguments:
being the permutation group of the four time points. Since the two-particle Green’s function depends on four discrete time arguments, it is not feasible to store in computer memory. We therefore evaluate Eq. (35) on the fly every time it is needed in the perturbation expansion. This can be done most efficiently if one precalculates the time-dependent creation and annihilation operators and only performs the matrix product and the subsequent trace at each time combination. Since this procedure is based on standard matrix multiplication, the computational effort of a single evaluation of Eq. (35) scales as , where is the dimension of the largest symmetry block in the Hamiltonian. An alternative way to calculate the ED quantities is to calculate the time dependence of the eigenstates from the Schrödinger equation and storing the operator matrices in memory as a function of time Eckstein and Werner (2010). The diagrams can then be evaluated using higher-order quadrature rules. Depending on the number of time points it is possible to efficiently treat reference systems up to a total system size of at least four sites in total.
Iii Application of the method
iii.1 Simple test
As a first test we calculate the time evolution of an exactly solvable model. Figure 3 shows the model under consideration.
The full system consists of an interacting site coupled to one additional bath site. At time the system is prepared in such a way that both sites are half filled and the spin degeneracy on the interacting site is lifted via a small magnetic field. The time dependence of the full system consists of a sudden switch in the hopping amplitude to a non-zero value.
The reference system, which is used as a starting point for the perturbation expansion, is modelled in the same way, but the hopping is switched to a different value. At this point we have to mention that although the auxiliary Hamiltonian used for preparing the initial state is correlated, an extension of the contour to imaginary times is not necessary. Here has exactly the same form as Eq. (8) and can therefore be treated on the contour depicted in Fig. 1. The reason for this is that the correlated impurity is decoupled from the bath for times smaller than zero. Consequently the density matrix splits into a direct product of a bath and impurity part. Fig. 3 shows the time dependence of the occupation number on the interacting site. The black dotted curve is the exact time evolution of the full system, the blue curve is the time evolution of the reference system. The green and red data points show different expansion orders in the dual potential. The green points correspond to a dual theory without any diagram, the red curve to the solution including the first diagram. Since this is a finite system, the solution is periodic (the period is larger than the time interval shown in the plot). The period of the oscillations of the magnetization is determined by the hopping amplitude . As one can see, both approximations improve the solution of the reference system towards the solution of the full system. In particular, already the zero-order solution corrects the oscillation period of the reference system. The quality of the solution increases with the the order of approximation. The solution including the first diagram gives the best improvement. The superperturbation Keldysh theory is in quite good agreement with the exact result.
In order to eliminate the systematic discretization error in the time argument, simulations for different grid sizes have been performed and the limit to a continuous time variable has been done numerically by quadratic regression. The final result can be expanded into a Taylor series in around the continuous solution:
The three constants , and the solution for a continuous time, , have been calculated by quadratic regression. The procedure is depicted in Fig. 4 for the time point . Here the dependence on is almost linear, but the quadratic regression is necessary to properly resolve the effect of the small initial magnetic field. The differences in for different are shown in Fig. 4.
iii.2 Spin in fermionic bath
An important non-trivial test for the method is the dissipation of a single spin, into and an infinitely long fermionic chain. Although the solution of the reference system is necessarily periodic, for the superperturbation solution of the infinite system, this may not be so.
The system under consideration consists of an isolated single magnetic Hubbard site which is coupled to the infinite chain at time by a sudden switch in the hopping to the chain. Since the chain is infinite, the problem is a very simple example for a model coupled to an infinite fermionic reservoir. In the following we want to demonstrate that although the dual perturbation expansion starts from a Hamiltonian system, which obeys energy conservation, dissipation can be found already in the lowest order approximation, without any dual diagram. This zero order approximation is very similar to a time-dependent Hubbard-I expansion (in the case , for the equilibrium case see e.g. Ref.Krivenko et al. (2010)) or a time-dependent variational cluster approach (t-VCA) (when , for equilibrium case see e.g. Ref.Potthoff et al. (2003)) and has the following very simple formulation:
Here is the approximation for the Green’s function of the full system and the Green’s function of the reference system. For a chain the hybridization function can be defined in terms of a continued fraction:
where is the non-interacting Green’s function of a single fermionic site and a time-dependent hopping matrix. Both quantities have exactly the same form as shown in Example 1. The details of the model and results are shown in Figure 5. The reference system again consists of the same Hubbard impurity coupled to a single bath level.
At time zero a single spin is prepared on the impurity, then a hopping to the rest of the chain is switched on.
Figure 5 shows the time-development of the magnetization on the impurity for the two site reference system (blue curve), in zero order dual perturbation theory (red curve) and of an impurity with the same parameters and a 5 site chain (green curve). At small times the length of the chain does not matter for the electron propagation, so all three curves coincide. At later times clear deviations are visible: The time-development of the reference system is governed by fast oscillations, which are caused by the reflection of the electron at the ends of the chain. The time-development of the magnetization for the model with 5 sites in the chain also shows a revival peak caused by the reflection of the electron, but this peak is found at later times, because of the longer chain.
The result of the perturbation expansion shows no backscattering peak and oscillates around zero, which is a clear indication that the presented approach allows to treat open systems starting from a perturbation around a closed Hamiltonian system.
It is clear that the perturbation expansion gives the best results, if the reference system is chosen in an optimal way, which essentially means that one should minimize a predefined matrix norm of by varying the effective parameters of the reference system. In the example at hand the minimum the effective hopping parameter of the reference system was chosen the same as the hopping parameter of the full system. The possibility to vary the reference system in principle allows to use the present solver as a solver for time-dependent DMFT calculations.
iii.3 Comparison to CT-QMC
In the next step, we compare the dual approach to the recently developed CT-QMC method for non-equilibrium systems Werner et al. (2009). The underlying model is a single-level quantum dot with Coulomb interaction . The spin-degenerate level with orbital energy is coupled to two fermionic reservoirs (, ) by a hybridization . The model Hamiltonian is given by
where represents the isolated dot with dot operators :
Furthermore, and describe a free electron gas in the left and right reservoir
where are creation and annihilation operators for the conductance electrons with momentum and spin . The energy distribution of the electrons in the two reservoirs obeys the Fermi function , assuming each reservoirs to be in thermal equilibrium with corresponding electro-chemical potentials and . A bias voltage can be applied to the dot by shifting the potentials . However, in the following we consider only the case with . The dot is symmetrically coupled to the left and right leads by the spin-conserving tunneling Hamiltonian
The strength of the tunneling can be characterized by the constant , where is the density of states in the leads. The density of states is taken to be constant with a hard cutoff Werner et al. (2009) at . In what follows, the energy scales are defined by the level broadening .
We introduce the reference system for the dual approach resembling the original model but reducing the leads to a single electronic level
The energy of the left and right lead is chosen to be equal to the orbital energy of the dot. For the model we investigated, the best agreement with MC is achieved when both chemical potentials in the auxiliary system equal the value of the level energy . In contrast, the chosen level energy of the axillary system related to the level energy of the dot is of minor importance.
Fig. 6 shows the results of the dual approach in lowest order (lines) for different values of compared to data obtained by the diagrammatic Monte Carlo method (symbols). As the non-equilibrium Monte Carlo method suffers from the dynamical sign problem we restrict ourselves to the short-time dynamics of the system. The electron density was calculated for the specific parameters leading to the half filled case with and temperature . Initially, the dot is empty and the coupling between the dot and the leads is switched on at . The comparison between the dual approach and non-equilibrium Monte Carlo shows good agreement for the chosen parameter range. For very strong interaction and large times, however, higher order corrections of the dual perturbation series become significant.
In the presented work we have generalized the dual fermion method to treat strongly correlated systems out of equilibrium. We have formulated the dual perturbation theory on the Keldysh contour and described the implementation of lowest order dual-diagram for a quite broad class of physical non-equilibrium problems. We also presented the details of the numerical implementation within the exact diagonalization scheme. First tests show that the method works as well for finite as for open systems in a fermionic bath. The long timescale non-equilibrium propagation proved to be accessible within the dual fermion perturbation theory, which is an advantage compared to non-equilibrium CT-QMC based methods. We also considered a simpler version of the presented method (a zero-order approximation), the time-dependent variational cluster approximation, that demands much less numerical effort and seems to be a promising tool for more complex systems. The dual fermion non-equilibrium scheme allows to go beyond zero-order and include physically important diagrammatic series which can describe the long time scale propagation of quantum systems out of equilibrium.
We thank M. Eckstein, I. Krivenko, M. Balzer, M. Potthoff, Ph. Werner, and A. Cavalleri for valuable discussions. This work has been supported by the Cluster of Excellence Nanospintronics (LExI Hamburg) and DFG Grant(436113/938/0-R).
- S. Wall, D. Prabhakaran, A. Boothroyd, and A. Cavalleri, Phys. Rev. Lett. 103, 97402 (2009).
- R. Tobey, D. Prabhakaran, A. Boothroyd, and A. Cavalleri, Phys. Rev. Lett. 101, 197404 (2008).
- S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005).
- R. Landauer, Phil. Mag. 21, 863 (1970).
- H. Schoeller and G. Schön, Phys. Rev. B 50, 18436 (1994).
- N. Wingreen and Y. Meir, Phys. Rev. B 49, 11040 (1994).
- K. Thygesen, Phys. Rev. Lett. 100, 166804 (2008).
- A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
- P. Werner, T. Oka, M. Eckstein, and A. Millis, Phys. Rev. B 81, 35108 (2010).
- A. Rubtsov, M. Katsnelson, and A. Lichtenstein, Phys. Rev. B 77, 33101 (2008).
- H. Hafermann, G. Li, A. N. Rubtsov, M. I. Katsnelson, A. I. Lichtenstein, and H. Monien, Phys. Rev. Lett. 102, 206401 (2009a).
- H. Hafermann, C. Jung, S. Brener, M. Katsnelson, A. Rubtsov, and A. Lichtenstein, EPL 85, 27007 (2009b).
- L. V. Keldysh, Sov. Phys. JETP 47, 1515 (1964).
- L. V. Keldysh, Sov. Phys. JETP 20, 1018 (1965).
- M. Wagner, Phys. Rev. B 44, 6104 (1991).
- A. Kamenev and A. Levchenko, Adv. Phys. 58, 197 (2009).
- J. W. Negele and H. Orland, Quantum many-particle systems (Perseus Books, Reading, MA, 1998).
- J. K. Freericks, Phys. Rev. B 77, 075109 (2008).
- M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010).
- I. S. Krivenko, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, JETP Let. 91, 339 (2010).
- M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. 91, 206402 (2003).
- P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 035320 (2009).