Hubbard nanoclusters far from equilibrium

Hubbard nanoclusters far from equilibrium

Sebastian Hermanns    Niclas Schlünzen    Michael Bonitz Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, Leibnizstraße 15, D-24098 Kiel, Germany
September 3, 2019

The Hubbard model is a prototype for strongly correlated many-particle 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 non-perturbative 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 HF-GKBA cures some artifacts of prior two-time 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.


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 out-of-equilibrium dynamics are of great current interest in solid-state, atomic and molecular physics, in nanoelectronics, quantum transport etc. In all these fields, the availability of intense and coherent radiation, combined with ultra-short 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 multi-electron atoms and molecules, e.g. [becker12, , schuette_prl12, ] and references therein. Another example are the many-body 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 many-body systems in nonequilibrium, including the optical excitation of electron-hole 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 carrier-phonon 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 one-dimensional (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 HF-GKBA 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 (HF-GKBA) retains the conserving properties of the original selfenergy. The constraints on the propagators are reformulated giving rise to a generalization and relaxation of the well-known energy conservation criteria of Baym and Kadanoff baym61 ; book_kadanoffbaym_qsm .

We, furthermore, report excellent numerical behavior (long-time stability of the propagation) of the HF-GKBA and confirm the advantageous scaling with the propagation time (, compared to in full two-time NEGF simulations [hermanns12_pysscripta, ]). Comparing to exact diagonalization results, we conclude that the HF-GKBA with second Born selfenergies is very accurate for weak coupling, for a time duration substantially larger than that of time-dependent Hartree–Fock and that this time range increases with the particle number. Finally, we use the HF-GKBA to study the short-time 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 single-time 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 many-particle 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 one-particle Green function, defined as (we set )


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.

Figure 1: Schwinger/Keldysh contour . The -branch starts from running in positive direction up to , the -branch starts from and leads in negative direction to . Note that on the contour the marked time is later than , although the corresponding times on the real time axis satisfy .

denotes a one-particle annihilation (creation) operator in an arbitrary one-particle 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 time-dependent Hamiltonian is given by


In second quantization, this expression attains the form


where denotes the one-particle Hamiltonian including an external potential, and is an arbitrary (possibly time-dependent) pair-interaction 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,


where , and is a delta function on the contour . The Martin–Schwinger hierarchy is equivalent to the exact many-body problem, describing the coupling of the evolution of the one-particle Green function to the two-particle Green function,


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 self-energy which is a functional of the single-particle Green function. This self-energy 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 :


which are the Keldysh–Kadanoff–Baym equations (KBE). The KBE are—in principle—exact equations of motion of the many-body 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 many-body 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 long-time behavior. Among the higher order selfenergies, we mention the second(-order) Born (2B), GW or T-matrix 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 2B-exchange 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.

Figure 2: (a) Set of nonequilibrium Feynman diagrams used in the present work, from left: Hartree, Fock and second order Born (direct and exchange) selfenergy. Wiggly lines denote the interaction potential, full lines are two-time Green functions. (b) When the HF–GKBA is applied, all full lines are replaced via Hartree–Fock Green functions Eqs. (17, 18), (thin lines). (c) Non-vanishing diagrams used in the HF-GKBA calculations for the Hubbard clusters. The exchange diagrams (gray lines) do not appear.

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 two-time plane by computing the NEGF as a function of both time arguments. Due to the time-memory 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 long-time 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 one-particle 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 real-time 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:


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 Langreth-Wilkins rules, e.g. [LWR, , balzer_lnp13, ]. For the correlation functions, we have


where the components “R/A” and “” of are defined analogously to those of G. The propagators satisfy


Lipavský et al. [lipavsky86, ] have shown that the KBE for are equivalent to an integral equation ():


whereas for :


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 two-time 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:


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 off-diagonal Green functions (i.e. for arguments ) from single-time 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 two-time 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 time-ordering operator):


where denotes the Hartree–Fock (mean–field) Hamiltonian,


that contains the time-dependent 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 single-particle 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 T-matrix 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


where denotes the ideal Green function of the uncorrelated and field-free 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 split-off the correlation selfenergy and introduce the HF-Green function according to


We defined as the correlation selfenergy in which all two-time 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 two-time 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 two-time 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 two-time 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 many-body systems. This is, in particular, relevant for the many-body 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:


the single-particle Green function obeys simultaneously the KBE and its adjoint, Eqs. (4) and (5), and


the two-particle Green function is symmetric with respect to both particles, i.e.

These conditions easily allow one to select conserving approximations for the two-particle 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, ],


where, in the end, is set, and we introduced the correlation part of the two-particle 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 Bethe-Salpeter equation. In what follows, it will be sufficient to consider the screened ladder approximation (SCA), Ref. [bornath99, ; bonitz_jpcs13, ]


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 T-matrix 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 HF-approximation for the propagators (18) is made, and our numerical results for the HF-GKBA 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 non-Markovian equation for the time-diagonal element of the Green function, i.e. for the one-particle density matrix. Such an equation can, of course, be derived independently from a one-time 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 long-time behavior, can be used to analyze the properties of the single-time 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


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 (BBGKY-hierarchy) follow from the von Neumann equation for by taking the partial trace,

where the two-particle 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 two-time 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 three-particle density operator, see also Ref. [hermanns_jpcs13, ]. Since the density operator approach does not involve two-time Green functions, full agreement with NEGF theory requires, in addition, the time-diagonal 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 HF-GKBA 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


with all other terms left unchanged compared to the HF-GKBA.

It is this renormalization of the two-particle 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 one-particle and two-particle density operators,


We now compute the time derivative of kinetic energy using the first hierarchy equation,


The time derivative of the potential energy follows similarly,


Finally, the time derivative of the interaction energy is transformed using the second hierarchy equation with the replacement (30)


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,


For a conserving approximation, the right hand side has to be equal to zero. The first term vanishes if the three-particle 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 two-particle Green function, cf. Sec. III.3. Here, in the case of single-time operators, it appears on the three-particle 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 pair-correlation operator, details can be found in Refs. [book_bonitz_qkt, ] and [hermanns_jpcs13, ],


where , is the two-particle Hartree–Fock hamiltonian and is the three-particle (anti-)symmetrization operator. Obviously, this set of equations assures symmetry of in the particle indices. This approximation is equivalent to the HF-GKBA in second Born approximation, i.e. no renormalization of the two-particle hamiltonian (30) is performed.

As noted above, for the HF-GKBA, the second term on the r.h.s. in Eq. (36) is absent and energy conservation is confirmed for any conserving approximation of two-time 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 anti-hermitean 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 single-particle propagator is renormalized by a single-particle 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 HF-GKBA 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 real-time 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 two-particle Green function is eliminated in favor of the selfenergy Keldysh matrix with the same components and and the same link between them


In our argument in Sec. III.3 we used the condition that fulfills its pair of KBE. But what about ? In standard two-time 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 HF-GKBA uses just , whereas the approximation with the exponentially damped propagators is associated with where is time-independent. Finally, we can establish a connection with the corresponding approximation for the two-particle 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 non-symmetric function . This explains the preservation of energy conservation for the HF-GKBA 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 time-dynamics of a system described by the single-particle Green functions and connected via a functional relation , with


obeying simultaneously the real-time KBE involving the selfenergy , Eq. (39), and its adjoint,


obeying simultaneously the real-time KBE involving the selfenergy , Eq. (40), and its adjoint,


and the two-particle Green functions associated with and being symmetric with respect to both particles, i.e.

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 single-particle 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 three-fold:

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

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

  3. 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 half-filling (particle number ) with hopping amplitude and on-site interaction . The initial Hamiltonian, for times , reads


where and label the discrete sites, and indicates nearest-neighbor sites. Further, denotes the density operator. The last term in Eq. (41) incorporates an external time-dependent 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 non-interacting 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 Fermi-like switching function


with a switching duration and a switching constant , which yields sufficiently converged results.

iv.2 Testing the GKBA for small clusters

Figure 3: (Color online) Top. Time evolution of the density on site for a 2-site Hubbard model at half filling and within HF–GKBA (full red line), compared to the exact result (full black line with dots) and a two-time second Born calculation (green dashed line). The HF–GKBA does not exhibit artificial damping, in contrast to the two-time result. Bottom. The same for a cluster of sites.

We start by considering small one-dimensional 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 two-time KBE simulations when the system was driven far out of equilibrium [puigvonfriesen09, ]. In that reference, at time , a two-site 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 two-time simulation in second Born approximation, while the exact result exhibits a non-decaying dynamics. The HF–GKBA, interestingly, does not exhibit the damping of the two-time result. Since the many-body 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 T-matrix 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.

Figure 4: (Color online) Top. Density evolution on the left-most site for an 8-site Hubbard model at half filling and coupling strength . Initially all particles were on the leftmost four sites. The HF–GKBA result (full red line) is compared to the exact one (full black line with dots) and to Hartree–Fock (dashed blue line with crosses). Bottom. Same as above, but for a reduced time interval. Instead of TDHF, a two-time second Born calculation (green dashed line) is included.

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 8-site Hubbard chain at weak coupling, , starting from a strong nonequilibrium situation where all particles are confined to the four left-most 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 HF-GKBA to full two-time 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 time-dependence 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 HF-calculation, 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.

Figure 5: (Color online) Energy spectrum for an 8-site Hubbard model at half filling and coupling strength . The HF–GKBA result (full red line) is compared to the exact one (full black line with dots) and to Hartree–Fock (dashed blue line with crosses).

After having tested the long-time behavior of the HF–GKBA in the linear response regime, let us now consider the long-time behavior in the case of a strong non-perturbative excitation. The results for a four-site 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 .

Figure 6: (Color online) Long-time dynamics of a 4-site Hubbard model at half filling and , following strong excitation, within HF–GKBA (red), compared to the exact result (black with dots). Initially all electrons occupy the two leftmost sites. The figures shows the time evolution of the occupation of the leftmost site.

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).

Figure 7: (Color online) Time evolution of the density on the leftmost site for an 16-site Hubbard model at half filling and within HF–GKBA (full red line), compared to the HF-result (blue dashed line with crosses). Initially all particles were at the 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 Short-time dynamics: correlation build up and relaxation of site occupations

Figure 8: Time evolution of the 8-site Hubbard model at half filling and . Initially, all particles were at the leftmost sites. Top figure: exact diagonalization. Bottom figure: HF–GKBA result. Upper panels of both figures show the dynamics of kinetic, Hartree–Fock and correlation energy, bottom panels of both figures show the occupation of the four leftmost sites.

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 2-4 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 .

Figure 9: Same as Fig. 8 but for . Top figure: Exact diagonalization calculation. Bottom figure: HF–GKBA

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 .

Figure 10: Same as Fig. 8 but for . Exact diagonalization calculation.

We now turn to larger particle numbers. Here, no exact d