Hubbard nanoclusters far from equilibrium
Abstract
The Hubbard model is a prototype for strongly correlated manyparticle systems, including electrons in condensed matter and molecules, as well as for fermions or bosons in optical lattices. While the equilibrium properties of these systems have been studied in detail, the nonequilibrium dynamics following a strong nonperturbative excitation only recently came into the focus of experiments and theory. It is of particular interest how the dynamics depend on the coupling strength and on the particle number and whether there exist universal features in the time evolution. Here, we present results for the dynamics of finite Hubbard clusters based on a selfconsistent nonequilibrium Green functions (NEGF) approach invoking the generalized Kadanoff–Baym ansatz (GKBA). We discuss the conserving properties of the GKBA with Hartree–Fock propagators in detail and present a generalized form of the energy conservation criterion of Baym and Kadanoff for NEGF. Furthermore, we demonstrate that the HFGKBA cures some artifacts of prior twotime NEGF simulations. Besides, this approach substantially speeds up the numerical calculations and thus presents the capability to study comparatively large systems and to extend the analysis to long times allowing for an accurate computation of the excitation spectrum via time propagation. Our data obtained within the second Born approximation compares favorably with exact diagonalization results (available for up to 13 particles) and are expected to have predictive capability for substantially larger systems in the weak coupling limit.
pacs:
I Introduction
Strongly correlated quantum systems and materials, e.g. [pavarini11, ], are of rapidly growing relevance in many fields of physics and chemistry. Especially the outofequilibrium dynamics are of great current interest in solidstate, atomic and molecular physics, in nanoelectronics, quantum transport etc. In all these fields, the availability of intense and coherent radiation, combined with ultrashort laser pulses, has triggered many key experiments that allow one to investigate matter under extreme nonequilibrium conditions where strong correlations and nonlinear effects occur simultaneously [balzer_jpcs13, ]. Examples are the photoionization of multielectron atoms and molecules, e.g. [becker12, , schuette_prl12, ] and references therein. Another example are the manybody dynamics of particles in lattice systems in condensed matter, e.g. [kehrein_njp10, , eckstein_prl09, , eckstein_prb10, ] and optical lattices [bloch08, ] following a rapid quench of the interaction strength.
From the theoretical point of view, such systems pose particular challenges since quantum, spin and strong correlation effects have to be treated selfconsistently under situations far from the ground state or from thermodynamic equilibrium. Here, remarkable progress has been achieved recently using ab initio methods such as exact diagonalization, density matrix renormalization group approaches, nonequilibrium dynamical mean field theory [eckstein_prl09, ], iterative path integral techniques [thorwart_13, ] and others. The common problem of these approaches is an exponential scaling with the particle number and restrictions with respect to the duration of the time propagation.
For this reason, alternative approaches that are based on statistical methods are of high interest. This includes density operator methods, e.g. [book_bonitz_qkt, ; akbari_12, ; hermanns_jpcs13, ], nonequilibrium Green function (NEGF) techniques and a recently developed stochastic mean field approach [ayik09, ; lacroix14, ]. Here, we focus on the NEGF method which, during the past 15 years, has been successfully applied to a variety of manybody systems in nonequilibrium, including the optical excitation of electronhole plasmas in semiconductors [kwong98, , kwong00, ], nuclear collisions [rios11, ], dynamics of laser plasmas [bonitz_cpp99, , haberland_pre01, ] and the problem of baryogenesis in cosmology [garny11, ]. More recently, NEGF methods have also been used to describe finite spatially inhomogeneous systems, including the carrier dynamics and carrierphonon interaction in quantum dots and quantum wells [gartner06, ; lorke06, ; bonitz_prb07, ; balzer_prb09, ], molecular transport in contact with leads, e.g. [uimonen11, ; khosravi12, ; stefanucci_gkba, ], or small atoms, e.g. [dahlen07, , balzer_pra10, , balzer_pra10_2, ]. For a recent overview on NEGF applications to inhomogeneous systems, see Refs. [balzer_lnp13, , stefanucci_book_13, ].
Applications of NEGF methods to small Hubbard clusters have been presented not long ago [puigvonfriesen09, ; puigvonfriesen10, ] and showed the great potential of this method. The physical features that could be explored include the relaxation dynamics, the excitation spectrum and, in particular, the relevance of double excitations [balzer12_epl, ; sakkinen12, ]. At the same time, NEGF simulations exhibited fundamental problems: The first is of conceptual nature and is related to unphysical damping effects in case of strong excitation which are absent in exact calculations [puigvonfriesen09, ]. The second is the computational difficulty due to the strong increase of CPU time and memory demand with increased propagation time which limits the spectral resolution and the duration of nonequilibrium calculations. We have recently presented an idea how to overcome the first and substantially weaken the second problem: invoking the generalized Kadanoff–Baym ansatz (GKBA) of Lipavský et al. [lipavsky86, ]. This concept was tested on the level of second order Born selfenergies for the example of a onedimensional (1D) Hubbard cluster containing just two sites and two electrons because here comparisons with available exact diagonalization methods are easily possible, cf. [hermanns_jpcs13, , hermanns12_pysscripta, , bonitz_cpp13, ].
Based on these encouraging first results, in this paper, we present a systematic analysis of the HFGKBA approach in application to Hubbard nanoclusters. We discuss fundamental questions such as the issue of total energy conservation and we extend our previous results to larger systems. This issue has regained importance since the idea to use the GKBA in NEGF calculations has recently been taken up for charge dynamics in molecular junctions [stefanucci_gkba, ] and the dynamics of localization in 1D Hubbard chains [lev_14, ]. While Bar Lev et al. [lev_14, ] used Hartree–Fock propagators Stefanucci et al. [stefanucci_gkba, ] used the GKBA in combination with Hartree–Fock as well as damped propagators. Here, we demonstrate that only the use of Hartree–Fock propagators (HFGKBA) retains the conserving properties of the original selfenergy. The constraints on the propagators are reformulated giving rise to a generalization and relaxation of the wellknown energy conservation criteria of Baym and Kadanoff baym61 ; book_kadanoffbaym_qsm .
We, furthermore, report excellent numerical behavior (longtime stability of the propagation) of the HFGKBA and confirm the advantageous scaling with the propagation time ( , compared to in full twotime NEGF simulations [hermanns12_pysscripta, ]). Comparing to exact diagonalization results, we conclude that the HFGKBA with second Born selfenergies is very accurate for weak coupling, for a time duration substantially larger than that of timedependent Hartree–Fock and that this time range increases with the particle number. Finally, we use the HFGKBA to study the shorttime dynamics of Hubbard clusters and present applications to larger systems.
This paper is organized as follows: In Sec. II, we recall the basics of the nonequilibrium Green functions approach. The transition to a singletime description with the help of the GKBA is discussed in Sec. III, where we also discuss its conserving character, in dependence on the choice of the propagators. Our numerical results for finite Hubbard clusters are presented in Sec. IV.
Ii Nonequilibrium Green functions
To describe correlation effects and excitations in quantum manyparticle systems, the NEGF approach has proven very successful, as it allows for a systematical inclusion of correlations by diagrammatic expansions. In contrast to density matrix based schemes, the Green function method additionally offers direct access to the spectral information as well as particle removal and addition energies. The main quantity is the oneparticle Green function, defined as (we set )
(1) 
where the brackets denote thermodynamic averaging and is the time ordering operator on the Schwinger–Keldysh round trip contour [KeldyshContour, ] , on which the times and are defined, see Fig. 1.
denotes a oneparticle annihilation (creation) operator in an arbitrary oneparticle basis in second quantization. To simplify the notation we suppress the orbital index “i” regarding as vectors . Correspondingly, the Green function (1) is understood as a matrix with components .
For a quantum system of particles in nonequilibrium, the generic timedependent Hamiltonian is given by
(2) 
In second quantization, this expression attains the form
(3) 
where denotes the oneparticle Hamiltonian including an external potential, and is an arbitrary (possibly timedependent) pairinteraction potential [here and are understood as matrices and , respectively].
The equations of motion for are the first equation of the Martin–Schwinger hierarchy [martin_theory_1959, ] and its adjoint,
(4)  
(5)  
where , and is a delta function on the contour . The Martin–Schwinger hierarchy is equivalent to the exact manybody problem, describing the coupling of the evolution of the oneparticle Green function to the twoparticle Green function,
(6) 
which itself is coupled to the three–particle Green function by a similar equation (the Bethe–Salpeter equation) and so on.
To solve this hierarchy and to make it numerically tractable, a formal decoupling is performed by introducing the selfenergy which is a functional of the singleparticle Green function. This selfenergy can be found from a diagrammatic expansion in terms of Feynman diagrams, where only some classes of diagrams are chosen according to the properties of the examined system. With this, the equations of motion (4) and (5) become formally closed equations for :
(7)  
(8) 
which are the Keldysh–Kadanoff–Baym equations (KBE). The KBE are—in principle—exact equations of motion of the manybody system would the selfenergy be exactly known. This is the case only for a limited number of models. In general, therefore, one has to resort to manybody approximations for the selfenergy. Baym and Kadanoff have shown how to select approximations that obey the conservation laws [baym61, ; book_kadanoffbaym_qsm, ], which we will discuss in Sec. III. The simplest approximation is the Hartree–Fock (HF) approximation where correlations are neglected entirely. It is commonly expected that this is a reasonable approximation for weak coupling. Nevertheless, we will see that, even for small coupling strength, in some nonequilibrium situations correlation effects may play a crucial role, in particular, for the longtime behavior. Among the higher order selfenergies, we mention the second(order) Born (2B), GW or Tmatrix approximations [book_bonitz_qkt, ]. In this paper, we will focus on the second order Born approximation shown diagramatically in Fig. 2, as it allows for long time simulations. Note that, due to the particular nature of the interaction in the Hubbard model, the Fock and the 2Bexchange diagrams (gray, second and fourth) vanish. For the treatment of Hubbard nanoclusters in higher order approximations, we refer to Ref. [puigvonfriesen10, ]. The remarkable property of the NEGF theory is that, via utilization of the roundtrip contour [KeldyshContour, ], all approximations known from ground state theory and thermodynamic equilibrium situations remain fully valid in nonequilibrium, including slow and rapid processes as well as weak and strong excitation.
The direct numerical solution of the KBE (7,8) is now routine, e.g. [book_bonitz_qkt, , book_bonitz_semkat, , balzer_lnp13, , stefanucci_book_13, ] and references therein. After preparing a correlated initial state, e.g. [semkat_jmp00, , stefanucci13, ], the system is propagated in the twotime plane by computing the NEGF as a function of both time arguments. Due to the timememory structure of the collision integral in Eqs. (7, 8), the NEGF at all times and for all values of the orbital (site and spin) indices have to be stored in memory [fedvr, ]. Here, substantial advances could be recently achieved via a sophisticated program structure and parallelization [balzer_pra10, ; balzer_pra10_2, ]. Nevertheless, the computational requirements for the KBE solutions exhibit an unfavorable cubic scaling with time [hermanns12_pysscripta, ]. Clearly, this limits the duration of propagation in nonequilibrium as well as the accuracy and resolution of the computed energy spectra that are obtained from a Fourier transform (time integral over the whole simulation). To overcome these limitations and make the longtime calculations feasible, we will apply the generalized Kadanoff–Baym ansatz (GKBA) which leads to a quadratic scaling with time and also has a number of other attractive features. This ansatz is discussed in detail in the next section.
Iii The generalized Kadanoff–Baym ansatz
The oneparticle Green functions appearing in Eqs. (4) and (5) depend on two times both of which can be located on either the upper or lower branch of the contour , cf. Figure 1. We note that we do not use a contour with a third (vertical imaginary) branch to produce a correlated initial state, e.g. [stefanucci13, ], [balzer_lnp13, ]. Instead, we will prepare this state by an initial realtime propagation during which the interaction is turned on adiabatically, e.g. [rios11, ; hermanns12_pysscripta, ]. Thus on represents a matrix, i.e. 4 functions with two conventional real time arguments . Two of these functions are independent. It is common to use the following definitions for the correlation () and retarded (R) and advanced (A) functions:
(9)  
(10)  
(11) 
To make the relations to the field operators clear, we temporarily restored the orbital indices . In the following, these indices will be suppressed again, i.e. and have to be understood as matrices and . The equations of motion for the correlation functions and the propagators follow directly from the KBE (7, 8) on the contour , applying the LangrethWilkins rules, e.g. [LWR, , balzer_lnp13, ]. For the correlation functions, we have
(12)  
(13)  
where the components “R/A” and “” of are defined analogously to those of G. The propagators satisfy
(14)  
Lipavský et al. [lipavsky86, ] have shown that the KBE for are equivalent to an integral equation ():
(15)  
whereas for :
(16)  
Note, that by exchanging () and replacing the density matrix by in equations (15) and (16), the analogous expression for is easily obtained. Details of the derivation can be found in Ref. [hermanns12_pysscripta, ].
iii.1 Gkba
While expressions (15) and (16) are still exact (within the chosen approximation for ) they contain the unknown twotime function also under the integral on the r.h.s. Therefore, one can attempt to solve the integral equation up to second order, approximating under the integral just by the first term:
(17) 
This is the generalized Kadanoff–Baym ansatz (GKBA) of Lipavský, pika and Velický which is exact on the time diagonal, . The importance of this equation lies in the fact that it provides a means for the reconstruction—though approximately—of the offdiagonal Green functions (i.e. for arguments ) from singletime quantities such as the density matrix . Note that the argument of is not the mean of the two times appearing on the left but always the earlier of the two times. This means the GKBA retains the retardation structure (memory) of the collision integrals and thus obeys causality, which turns out to be crucial for the conservation properties, see Sec. III.3.
Note that Eq. (17) indicates that the GKBA is only formally closed in terms of , since it still involves twotime quantities—the retarded (advanced) propagators . These functions obey equations of motion of similar complexity as , cf. Eq. (14). To make further progress, one can use approximate propagators that are not computed selfconsistently with .
iii.2 Hartree–Fock GKBA
In this paper, we approximate the propagators by Hartree–Fock propagators, that are obtained from the solution of the KBE (14) with the replacement , with the solution ( is the causal timeordering operator):
(18) 
where denotes the Hartree–Fock (mean–field) Hamiltonian,
(19)  
(20) 
that contains the timedependent external field [via ] and interaction effects via the Hartree–Fock (mean field) selfenergy that selfconsistently involves the time–dependent density matrix. For later purposes, we also introduced the definition for the sum of Hartree–Fock selfenergy and external field, with being just the stationary singleparticle energy (kinetic plus potential energy).
Approximation (18), together with the GKBA, Eq. (17), will be called Hartree–Fock GKBA (HF–GKBA). One motivation of this choice is that, in the special case that the system is treated within the HF approximation (neglecting the correlation contribution, ), the HF–GKBA is exact, i.e. Eqs. (17, 18) provide the exact solution for which is readily confirmed by direct solution of the KBE (12, 13) in Hartree–Fock. Note that HF–GKBA can be used for an arbitrary correlation selfenergy . The example of the Tmatrix selfenergy will be briefly discussed in Sec. V.
Let us now have a closer look at the HF–GKBA for the case that correlations are taken into account and discuss its consequences. The Dyson equation for the full Green function on the contour can be written as
(21) 
where denotes the ideal Green function of the uncorrelated and fieldfree system. As discussed in Ref. [bonitz_cpp99_2, ], this equation can be decomposed into several coupled integral equations that allow to construct the full Green function in steps. Here, we choose to splitoff the correlation selfenergy and introduce the HFGreen function according to
(22)  
(23)  
(24) 
We defined as the correlation selfenergy in which all twotime functions are reconstructed according to the HF–GKBA, Eqs. (17, 18), thus, . The deviation of the full selfenergy from this approximation has been denoted and contains terms with one, two and three full propagators and , respectively. Thus, the Green function computed within the HF–GKBA is given by Eq. (23), containing renormalizations by the Hartree–Fock selfenergy and the second Born selfenergy with propagators on the Hartree–Fock level. In contrast, the full Green function that is computed by a full twotime calculation (i.e. without the GKBA) has undergone another renormalization given by Eq. (24). In addition to , the function contains contributions from to all orders, so the difference of the two is given by the infinite series
where each subsequent term adds contributions with up to three full (renormalized by correlations) propagators .
The properties of the GKBA have been studied for macroscopic spatially homogeneous systems [bonitz_jpcm96, ]. There, it was found that the GKBA retains the conservation laws of the original twotime approximation for the selfenergy [book_bonitz_qkt, ; bonitz_pla96, ]. This will be considered in more detail in Sec. III.3. As to the accuracy of the GKBA, it was found that this ansatz is a very good approximation to the full twotime solution if the exact propagators are being used, indicating that the additional integral contributions in Eqs. (15) and (16) are often of minor importance. In the case of approximate propagators, good results were obtained with ideal as well as Hartree–Fock propagators [kwong98, ]. In contrast, the use of damped propagators that include imaginary selfenergy contributions, violates total energy conservation and leads to an overall worse performance [bonitz_epjb99, ]. The use of the HF–GKBA for finite Hubbard clusters [hermanns12_pysscripta, ; bonitz_cpp13, ] confirms these results. Details will be presented in Sec. IV.
iii.3 Conservation of total energy
The issue of conserving approximations is of central importance for the treatment of correlated manybody systems. This is, in particular, relevant for the manybody dynamics far from equilibrium. It is an attractive feature of Green functions theory that conserving approximations are straightforwardly selected. Baym and Kadanoff [baym61, ; book_kadanoffbaym_qsm, ] have formulated a simple criterion for a NEGF approximation to be conserving that consists of two conditions:
 A
 B

the twoparticle Green function is symmetric with respect to both particles, i.e.
.
These conditions easily allow one to select conserving approximations for the twoparticle Green function and the selfenergy.
We now show that, when a conserving approximation for is being used, the subsequent application of the HF–GKBA does not change the conservation properties. Application of the GKBA amounts to solving the KBE only along the time diagonal, . The corresponding equation of motion for is obtained by computing the difference of the KBE and its adjoint, Eqs. (4, 5), cf. Ref. [book_kadanoffbaym_qsm, ],
(26) 
where, in the end, is set, and we introduced the correlation part of the twoparticle Green function, [recall that contains the Hartree–Fock selfenergy, cf. Eq. (20)]. By construction, the solution fulfills condition A.
Consider now condition B. To this end we use the solution for the correlation part of that follows from the BetheSalpeter equation. In what follows, it will be sufficient to consider the screened ladder approximation (SCA), Ref. [bornath99, ; bonitz_jpcs13, ]
(27)  
This equation can be solved by iteration, starting by replacing under the integral. This first iteration corresponds to the dynamically screened second Born approximation (GW) which, obviously, is symmetric in the labels of particles 1 and 2, in agreement with condition B. If we now apply the HF–GKBA to each Green function under the integral, this symmetry is fully retained. Thus, we have shown that the HF–GKBA for the GW approximation is conserving. The same applies to the static limit when , i.e. the static second Born approximation is conserving as well when the HF–GKBA is applied. The same proof applies to the Tmatrix approximation and to the SCA. To show this, we only need to proceed further with the iterative procedure for the solution of Eq. (27), either with the static or dynamic potential. It is easy to realize that each term of the iteration series has the needed symmetry , resulting in the fulfillment of condition B.
Thus, we conclude that the application of the GKBA to an arbitrary conserving approximation of NEGF theory does not change the exchange symmetry, condition B. This symmetry is also retained when, in addition, the HFapproximation for the propagators (18) is made, and our numerical results for the HFGKBA fully confirm total energy conservation. There is, however, a serious problem with the previous argument. Let us consider damped propagators, i.e. replace . This approximation is known to violate energy conservation [bonitz_jpcm96, ; bonitz_epjb99, ], although it also clearly obeys the symmetry and, thus, fulfills conditions A and B. In order to understand the origin of the violation of energy conservation for the exponentially damped propagators we, therefore, now first consider a different approach that is based on density operator theory. We will then return to the conditions A and B in Sec. III.5 and resolve this contradiction.
iii.4 Density operator theory and the GKBA
With the application of the HF–GKBA to the KBE, the problem becomes a closed nonMarkovian equation for the timediagonal element of the Green function, i.e. for the oneparticle density matrix. Such an equation can, of course, be derived independently from a onetime theory of reduced density matrices. This was first shown in Ref. [book_bonitz_qkt, ], see also Refs. [bonitz_jpcm96, ; hermanns_jpcs13, ]. This equivalence is important to identify the HF–GKBA with standard approximations from density operators. At the same time, results from density operator theory, including conservation laws and longtime behavior, can be used to analyze the properties of the singletime solutions of the KBE. Finally, we note recent interest in density operator methods in the context of the relaxation dynamics of finite Hubbard clusters [akbari_12, ].
We, therefore, briefly recall the concept of the reduced nonequilibrium density operators
(28)  
where is the density operator of the full system which is normalized to unity and is the associated particle operator. The equations of motion for the (BBGKYhierarchy) follow from the von Neumann equation for by taking the partial trace,
where the twoparticle hamiltonian is , and the system (LABEL:eq:bbgky) has to be complemented by initial conditions for etc. The equations are coupled and form a hierarchy that eventually stops when the r.h.s. involves , in analogy to the twotime Martin–Schwinger hierarchy of NEGF, see Sec. II. As in the case of NEGF, the hierarchy is usually decoupled at a low level by replacing the exact by an approximation that is a functional of the lower order operators. Key approximations of NEGF theory, including the second Born approximation [bonitz_jpcm96, ], ladder approximation [kremp_ap97, ] or GW approximation [book_bonitz_qkt, ] are readily identified by proper choices for the threeparticle density operator, see also Ref. [hermanns_jpcs13, ]. Since the density operator approach does not involve twotime Green functions, full agreement with NEGF theory requires, in addition, the timediagonal limit as provided by the GKBA, and in fact the GKBA is directly recovered in the theory of reduced density operators, e.g. [book_bonitz_qkt, ]. For the present purpose of analyzing energy conservation, it is sufficient to note that the HFGKBA is directly recovered from the system (LABEL:eq:bbgky). On the other hand, the GKBA with correlated propagators containing correlation selfenergy contributions [such as exponential damping] leads to a modification of the second hierarchy equation by the replacement
(30)  
with all other terms left unchanged compared to the HFGKBA.
It is this renormalization of the twoparticle hamiltonian that destroys the conservation of total energy in the GKBA with correlated propagators. To show this, we recall the derivation of conserving approximations in density operator theory [dufty_boercker, ; book_bonitz_qkt, ], for related issues of density operator theory, we refer to Refs. [dufty_boercker89, ; dufty_97, ]. We begin with expressing the mean kinetic, potential and interaction energy via the oneparticle and twoparticle density operators,
(31)  
(32) 
We now compute the time derivative of kinetic energy using the first hierarchy equation,
(33) 
The time derivative of the potential energy follows similarly,
(34)  
Finally, the time derivative of the interaction energy is transformed using the second hierarchy equation with the replacement (30)
(35)  
Collecting the results (33, 34, 35), we obtain, for the time derivative of the total energy minus the power introduced into the system by the external potential,
(36)  
For a conserving approximation, the right hand side has to be equal to zero. The first term vanishes if the threeparticle density operator is symmetric with respect to the particle indices, , at all times. This is an obvious and trivial condition [it is similar to condition B of NEGF theory for the twoparticle Green function, cf. Sec. III.3. Here, in the case of singletime operators, it appears on the threeparticle level], and is fulfilled also for the exact solution.
To verify this symmetry for the second Born approximation, we first rewrite in terms of correlation operators (cluster expansion) and give the corresponding expression for the paircorrelation operator, details can be found in Refs. [book_bonitz_qkt, ] and [hermanns_jpcs13, ],
(37) 
where , is the twoparticle Hartree–Fock hamiltonian and is the threeparticle (anti)symmetrization operator. Obviously, this set of equations assures symmetry of in the particle indices. This approximation is equivalent to the HFGKBA in second Born approximation, i.e. no renormalization of the twoparticle hamiltonian (30) is performed.
As noted above, for the HFGKBA, the second term on the r.h.s. in Eq. (36) is absent and energy conservation is confirmed for any conserving approximation of twotime NEGF theory. In contrast, for the GKBA with propagators that contain correlation selfenergies, , even for a conserving approximation (with a properly symmetric ) energy conservation is destroyed. In fact, just the antihermitean part of , which governs the damping behavior of the propagators, contributes to the second term on the r.h.s. This is particularly evident in the Born approximation where , and each singleparticle propagator is renormalized by a singleparticle correlation selfenergy contribution .
iii.5 Energy conservation of the GKBA revisited. Relaxing the conservation criteria of Baym and Kadanoff
After having confirmed energy conservation for the HFGKBA and violation thereof for damped propagators within an independent density operator approach, it remains to establish how this result can be obtained within NEGF theory. In particular, the question arises whether conditions A and B of Baym and Kadanoff are indeed sufficient and necessary for energy conservation.
We first find out why the symmetry of alone (condition B) does not imply energy conservation of the GKBA for arbitrary choice of the propagators as it appeared in Sec. III.3. Let us go back to condition A. The strict meaning of the statement (implied, by Baym and Kadanoff) that the Green function obeys simultaneously the KBE and its adjoint, Eqs. (4 ,5) is that each of the realtime Keldysh components simultaneously fulfills these equations. Since only two of these functions are independent it is sufficient to consider and which obey the pairs of equations (12) and (14), respectively. In these equations, the twoparticle Green function is eliminated in favor of the selfenergy Keldysh matrix with the same components and and the same link between them
(38) 
In our argument in Sec. III.3 we used the condition that fulfills its pair of KBE. But what about ? In standard twotime NEGF theory, of course, also fulfills its pair. But this is no longer the case when we apply the GKBA. In this approximation, the standard relation between and is altered and replaced by Eq. (17). This means, do not follow from the result for and but obey an independent equation. In other words, the equation of motion for , if written again in the standard form (14), will contain a selfenergy that is independent of and not given by Eq. (38). In fact, the GKBA just uses this independence in favor of a simpler choice for to simplify the calculations. For example, the HFGKBA uses just , whereas the approximation with the exponentially damped propagators is associated with where is timeindependent. Finally, we can establish a connection with the corresponding approximation for the twoparticle Green function from the standard mapping . While the Hartree–Fock selfenergy corresponds to which is obviously symmetric in the particle indices, the term is associated with a nonsymmetric function . This explains the preservation of energy conservation for the HFGKBA and its violation for the GKBA with exponentially damped propagators where both statements are independent of the original choice of provided it was conserving.
We may now use this result to revise the conditions A and B such that they apply to the GKBA. In fact we can relax the conditions of the theorem of Baym and Kadanoff (BKT) such that they cover not only the GKBA but a broader class of conserving approximations than envisaged originally in Refs. [baym61, ; book_kadanoffbaym_qsm, ].
Theorem: The timedynamics of a system described by the singleparticle Green functions and connected via a functional relation , with
 A1

obeying simultaneously the realtime KBE involving the selfenergy , Eq. (39), and its adjoint,
(39)  A2

obeying simultaneously the realtime KBE involving the selfenergy , Eq. (40), and its adjoint,
(40)  B

and the twoparticle Green functions associated with and being symmetric with respect to both particles, i.e.
and
,
is conserving.
Proof: Consider first Eq. (39). This equation is decoupled from . Thus, the symmetry condition B for guarantees, according to the BKT, that and, hence the dynamics of , are conserving. Consider now Eq. (40). It contains which, according to condition B, is conserving. Eq. (40), in addition, contains and . Since the dynamics of is conserving and the functional relation is a singleparticle relation having the same form for and , we conclude that the coupling to does not destroy the conservation properties of . Thus, the problem is reduced to the original BKT, and the dynamics of is conserving as well what proves the theorem.
Thus, instead of one selfenergy there are now two selfenergies that have to be conserving simultaneously. Obviously, this version of the theorem reduces to the BKT if the connection between and is given by the standard NEGF relation, Eq. (11) and, consequently, . If the relation is given by Eq. (17) we recover the GKBA where the choice of the retarded propagators is given by . Finally, there exists a whole set of new conserving approximations () where the functional relation (17) is replaced by a different one, .
Iv Results for finite Hubbard clusters
Let us now apply the NEGF formalism with the HF–GKBA to the dynamics of standard finite lattice systems described by the Hubbard model. The purpose of this section is threefold:

to test conceptual questions of the GKBA and the numerical performance. It has been reported before that twotime NEGF simulations for small Hubbard clusters exhibit unphysical damping behavior [friesen_10, ]. Furthermore, approximations based on density matrices have shown numerical instabilities in the longtime regime [akbari_12, ]. We will show in Sec. IV.2 that these problems do not occur with the HF–GKBA,

to extend the simulations to larger systems where no exact diagonalization data are available, cf. Sec. IV.2, and

to investigate the dynamical relaxation behavior of small Hubbard clusters at different coupling strengths, cf. Sec. IV.3.
iv.1 Hubbard model and second Born approximation
We consider a finite Hubbard model with sites at halffilling (particle number ) with hopping amplitude and onsite interaction . The initial Hamiltonian, for times , reads
(41) 
where and label the discrete sites, and indicates nearestneighbor sites. Further, denotes the density operator. The last term in Eq. (41) incorporates an external timedependent excitation that drives the system out of equilibrium. Several cases for the choice of the function will be specified below. In the following, we will use dimensionless parameters where energy (time) is measured in units of (the inverse hopping amplitude ) and the coupling strength and field amplitude will be given by and , respectively.
In our simulations, we start from an noninteracting initial state and adiabatically turn on the interactions to reach a fully correlated state. After this, the external excitation is turned on. Details on the procedure can be found in Refs. [rios11, ; hermanns_jpcs13, ; bonitz_cpp13, ]. We have verified in advance that the switching is slow enough to avoid any artifacts in the dynamics. In practice, we use a Fermilike switching function
(42) 
with a switching duration and a switching constant , which yields sufficiently converged results.
iv.2 Testing the GKBA for small clusters
We start by considering small onedimensional Hubbard clusters where exact diagonalization results are available. This allows for a rigorous test of the HF–GKBA results in its full time dependence and for a broad variety of excitation conditions.
Let us first investigate the problem of artificial damping of the dynamics of Hubbard clusters that was observed in full twotime KBE simulations when the system was driven far out of equilibrium [puigvonfriesen09, ]. In that reference, at time , a twosite Hubbard cluster at half filling is strongly perturbed by a rapid change of the energy of site “1”, which leads to the following choice for the last term in the Hamiltonian (41), cf. Ref. [balzer_jpcs13, ], , where . The results are shown in Fig. 3. One clearly sees the rapid damping of the density on the perturbed site, in the twotime simulation in second Born approximation, while the exact result exhibits a nondecaying dynamics. The HF–GKBA, interestingly, does not exhibit the damping of the twotime result. Since the manybody approximations are in both cases identical, the difference is solely due to the HF–GKBA and the broken selfconsistency, as was discussed in Sec. III.2, cf. Fig. 2. We note that we observe quantitative deviations from the exact result which, however, become substantially smaller when the particle number increases, cf. the figures below. The removal of the artificial damping is an important generic feature of the HF–GKBA and was confirmed in all our simulations. The same behavior is observed when the Tmatrix selfenergy is used, see Sec. V. This gives us confidence for applying this approximation to Hubbard clusters in nonequilibrium, in particular to larger systems where no exact diagonalization data are available.
Let us now turn to larger clusters and a more quantitative comparison with exact data. Figure 4 shows the dynamics of 8 electrons in an 8site Hubbard chain at weak coupling, , starting from a strong nonequilibrium situation where all particles are confined to the four leftmost sites by applying a strong confinement potential. At time , this potential is removed instantaneously and the dynamics are followed (this scenario was studied in Ref. [akbari_12, ]). One observes very strong oscillations of the site occupations as particles move to the right, towards the originally empty four sites. After about four periods, these oscillations become anharmonic and continue with a reduced amplitude. A Hartree–Fock calculation fails already after about two periods (compare the blue dashed curve to the exact result), indicating that the dynamics are strongly influenced by correlation effects. In contrast, our HF–GKBA approach yields very good agreement with the exact data for about 7 periods () after which deviations are increasing, but still most of the features are reproduced qualitatively. In particular, the dominant frequencies and peak positions are still captured. Comparing the HFGKBA to full twotime 2B results (green dashed line), we even see a better agreement with the exact solution throughout the whole simulation.
Next, we consider the energy spectrum for this system, again at weak coupling, . This is done by applying a weak very short external field pulse to site , , that excites all possible transitions. The HF–GKBA allows us to propagate the system for a long time, until , and to compute the timedependence of all observables. The energy spectrum is obtained via Fourier transform of the occupation of the perturbed site, , and the results resolve about seven orders of magnitude. There are two main energy regions. For frequencies below , all three approximations, including the Hartree–Fock simulation, show overall very good agreement. Only a few smaller features – the small side peaks of the peaks around and and the little peaks around are missing in the HFcalculation, indicating that these are correlation effects (most likely double exciations [balzer12_epl, ]). For frequencies , the picture changes completely. Evidently, HF misses all peaks. In contrast, the HF–GKBA performs impressively well, taking into account the very low height of the peaks in that range.
After having tested the longtime behavior of the HF–GKBA in the linear response regime, let us now consider the longtime behavior in the case of a strong nonperturbative excitation. The results for a foursite chain at half filling and are shown in Fig. 6. The first observation is that the simulations run stable for the entire duration of . This is confirmed by additional tests that are several times longer (not shown) indicating that previously reported stability problems [akbari_12, ] do not occur within the HF–GKBA. Further, the comparison with the exact result shows that the main frequencies are well reproduced, and even the phase of the oscillations is correct up to . Yet even at longer times the overall behavior is well captured, despite the dephasing, and deviations decrease again strongly, cf., e.g., the behavior around . At the same time, quantitative deviations are observed (amplitude of the oscillations), starting around .
These observations encourage us to extend our simulations also to larger systems where no exact results are available. In all cases we confirm that the HF–GKBA is completely stable. As an example, in Fig. 7, we show the dynamics of a 16 site system at half filling and with the same type of strong nonequilibrium initial conditions (all 16 particles occupy the 8 leftmost sites).
The behavior is similar to the dynamics of the previously studied analogous systems of four and eight particles, cf. Figs. 6 and 4, respectively. The site occupations undergo a rapid and violent evolution which is strongly influenced by correlation effects. Hartree–Fock simulations reproduce only the first periods of the main oscillation (). Based on our HF–GKBA results, we can deduce a number of trends when the particle number is increased: first, the main oscillation period increases proportional with and, second, the oscillations become increasingly nonlinear. Since the energy spectrum becomes much more complex when the system size increases it is presently not possible to relate this oscillation to a characteristic energy transition. To shed more light onto the physical processes involved in the dynamics and the dependence on the coupling strength and particle number, we will consider an increased set of quantities below, in Sec. IV.3.
Summarizing this first part of numerical results, we may conclude that the HF–GKBA is very well suited to study the dynamics of finite Hubbard clusters in the weak coupling regime, thereby (at least partially) overcoming problems of previous approaches. No unphysical damping and instabilities are observed. Since the present results are derived from selfenergies in second order Born approximation we have restricted the analysis to weak coupling, . For moderately larger values of the coupling parameter, still acceptable results can be obtained, as will be shown below for .
iv.3 Shorttime dynamics: correlation build up and relaxation of site occupations
In the following, we consider the dynamics in the same strong nonequilibrium situation that was discussed above, where all particles are initially placed on the leftmost sites (half filling). We now look at the behavior of additional observables. In particular, we follow the time dependence of the occupations of all initially occupied sites (the occupations of the rightmost sites follow by symmetry, e.g. , and so on).
The typical dynamics can be seen in the bottom panel of Fig. 8 displaying exact diagonalization results for the site chain at weak coupling, . Due to Pauli blocking, initially, only electrons from site can move to the right whereas electrons from site () can only follow when site () is being depopulated. This time delay in the depopulation is clearly visible. Interestingly, the rightmost site () is only depopulated half and sites and even less. Depopulation of the leftmost site () sets in last but proceeds to the lowest value of all sites (close to zero). This is easy to understand: while the population of sites 24 is increased again by newly incoming particles from the left, no such incoming flux exists at the boundary of the chain (site 1). After a very short time, , the rightmost site () is almost fully occupied, i.e., the electron wave has reached the right border after which it is being reflected. Subsequently, a strongly nonlinear oscillatory dynamics of the site occupations occurs that is damped until the occupations reach the stationary values of a homogeneous system (), around . This, however, is not a true stationary system, as our finite system has a reversible dynamics and, consequently, we clearly observe, at later times, a strong departure from the homogeneous configuration.
It is interesting to consider, besides the occupations, also the different contributions to the total energy [note that total energy is conserved to very high accuracy] of the system which are shown in the top panel of Fig. 8. In the initial state, the system has only mean field (Hartree–Fock) energy. Both, kinetic and correlation energy are exactly zero. When the occupied sites get depopulated, we observe a rapid increase of kinetic energy which is almost completely compensated by a loss of HF energy. Kinetic energy reaches its maximum (the particle current is largest) around after which it decreases again and continues to oscillate. It is interesting to note that this (first) maximum of kinetic energy is reached when the population of the leftmost is decreased to . This is observed in all simulations, independent of the coupling strength. Thus, this is a propagation effect resulting from the nonequilibrium inhomogeneous initial condition and constitutes the shortest time scale we observe in our simulations (“phase I”).
The second time scale that is apparent from the energy relaxation (“phase II”) is characterized by a formation of correlation energy and a decay and saturation of Hartree–Fock energy which is reached around . After this time, no qualitative changes of the energy dynamics are observed, aside from nonlinear oscillations of kinetic and correlation energy that occur with almost exactly opposite time derivatives. More characteristic for this phase III is the relaxation of the site occupations which terminates around , followed by a longer phase IV of occupation revivals, as mentioned above. In the following, we will analyze whether these four phases are visible also for other values of the coupling.
Consider now Fig. 9, where the dynamics for the same conditions are shown, but for the larger coupling strength . Again, we observe the rapid depopulation of the sites (phase I) and associated relaxation of kinetic and Hartree–Fock energy, with the same characteristic time . The main difference, compared to the case , is the more rapid saturation of HF energy and build up of correlation energy (phase II) terminating around . The relaxation of the occupations (phase III) lasts again until .
Finally, we consider the case , cf. Fig. 10. As before, phase I has a duration of . The relaxation of the occupations now takes slightly longer, until . The most striking difference to the previous cases, however, is in the dynamics of the correlation energy. Here, buildup of correlations is essentially over around whereas the Hartree–Fock energy saturates only around .
We now turn to larger particle numbers. Here, no exact d