Time-propagation of the Kadanoff-Baym equations for inhomogeneous systems
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J un l. Time-propagation of the Kadanoff-Baym equations for inhomogeneous systems ∗ Adrian Stan,
1, 2, 3
Nils Erik Dahlen, and Robert van Leeuwen
1, 2 Department of Physics, Nanoscience Center, FIN 40014, University of Jyv¨askyl¨a, Jyv¨askyl¨a, Finland European Theoretical Spectroscopy Facility (ETSF) Rijksuniversiteit Groningen, Zernike Institute for Advanced Materials,Nijenborgh 4, 9747AG Groningen, The Netherlands.
We have developed a time propagation scheme for the Kadanoff-Baym equations for generalinhomogeneous systems. These equations describe the time evolution of the nonequilibrium Greenfunction for interacting many-body systems in the presence of time-dependent external fields. Theexternal fields are treated nonperturbatively whereas the many-body interactions are incorporatedperturbatively using Φ-derivable self-energy approximations that guarantee the satisfaction of themacroscopic conservation laws of the system. These approximations are discussed in detail for thetime-dependent Hartree-Fock, the second Born and the GW approximation. PACS numbers: 31.15.xm, 31.15.-p
I. INTRODUCTION
The recent developments in the field of molecular elec-tronics have emphasized the need for further developmentof theoretical methods that allow for a systematic studyof dynamical processes like relaxation and decoherenceat the nanoscale. Understanding these processes is of ut-most importance for making progress in molecular elec-tronics, whose ultimate goal is to minimize the size andmaximize the speed of integrated devices [1]. To studythese phenomena, theoretical methods must allow for thepossibility to study the ultrafast transient dynamics [2, 3]up to the picosecond [4, 5] and femtosecond timescale,while including Coulomb interactions, without violatingbasic conservation laws such as the continuity equation[6]. A theoretical framework that incorporates thesefeatures is the nonequilibrium Green function approachbased on the real-time propagation of the Kadanoff-Baym (KB) equations [7, 8, 9, 10, 11, 12, 13, 14]. Thismethod allows for systematic inclusion of electron in-teractions while providing results in agreement with themacroscopic conservation laws of the system [6, 8]. In tworecent Letters [7, 11] we applied the KB equations to in-vestigate the short time dynamics of atoms and moleculesin time-dependent external fields, as well as the transportdynamics of double quantum dot devices. It is the aimof this paper to describe in detail the underlying methodthat was only briefly described in those Letters. This in-cludes both a description of the theory as well as the time-propagation algorithm. We further generalize the equi-librium method, described in two recent papers [15, 16],to the nonequilibrium domain. We also extend earlier ∗ See J. Chem. Phys. for the published version. work on the time-propagation method of the KB equa-tions for homogeneous systems [17, 18] to the case ofinhomogeneous systems. In the inhomogeneous case wecan not take advantage of Fourier transform techniquesanymore. The KB equations become time-dependent ma-trix equations instead, in which the matrices are indexedby basis function indices. The time-stepping algorithmhas to take into account the special double-time structureof the equations which are furthermore nonlinear, inho-mogeneous and non-Hermitian. Therefore, several stan-dard time-propagation methods can not be used. Our ap-proach is different from the one presented in Refs. [17, 18]by incorporating correlated initial states and the memorythereof, which is described in terms of Green functionswith mixed real and imaginary time arguments. To sim-plify the time-stepping procedure, we make use of severalsymmetry relations of the Green function.The paper is divided as follows: in section II we presentthe KB equations and their symmetry properties. In sec-tion III we discuss the conserving self-energy approxi-mations that we use, and in section IV we present thetime-propagation method that we developed for systemsdescribed within a general basis set representation. Fi-nally in section V we present a summary and conclusions.
II. THEORY
We consider a many-body system that is initially inequilibrium at a temperature T and with a chemical po-tential µ . At an initial time t the system is exposed to atime-dependent external field. This external field can, forinstance, be a bias voltage in a quantum transport case,or a laser pulse. The field forces the system out of equi-librium and we aim to describe the time-evolution of thisnonequilibrium state. In second quantization the time-dependent Hamiltonian of the system reads (throughoutthis paper we use atomic units ~ = m = e = 1)ˆ H ( t ) = Z d x ˆ ψ † ( x ) h ( x , t ) ˆ ψ ( x ) ++ 12 Z Z d x d x ˆ ψ † ( x ) ˆ ψ † ( x ) v ( r , r ) ˆ ψ ( x ) ˆ ψ ( x ) , (1)where x = ( r , σ ) denotes the space- and spin coordinates.The two-body interaction will, in general, be a Coulombicrepulsion of the form v ( r , r ) = 1 / | r − r | . The one-body part of the Hamiltonian is h ( x , t ) = − ∇ + w ( x , t ) − µ, (2)where w ( x , t ) is a time-dependent external potential.The chemical potential µ of the initial equilibrium sys-tem is absorbed in the one-body part of the Hamiltonian.The expectation value of an operator ˆ O , for a system ini-tially in thermodynamic equilibrium ( t < t ), is givenby h ˆ O i = Tr { ˆ ρ ˆ O } , (3)where ˆ ρ = e − β ˆ H / Tr e − β ˆ H is the statistical operator,ˆ H is the time-independent Hamiltonian that describesthe system before the time-dependent field is applied and β = 1 /k B T is the inverse temperature. The trace hererepresents a summation over a complete set of states inFock space [19]. After the time-dependent external fieldis switched on at time t , the expectation value is givenby h ˆ O ( t ) i = Tr { ˆ U ( t − iβ, t ) ˆ O H ( t ) } Tr { ˆ U ( t − iβ, t ) } , (4)where ˆ O H ( t ) = ˆ U ( t , t ) ˆ O ˆ U ( t, t ) is the opera-tor ˆ O in the Heisenberg picture and ˆ U ( t , t ) = T [exp( − i R t t dt ˆ H ( t ))], for t > t , is the time-orderedevolution operator of the system. We further wroteexp( − β ˆ H ) = ˆ U ( t − iβ, t ) as an evolution operator inimaginary time. If we read the time arguments in Eq.(4)from right to left we see that they follow a time-contouras displayed in Fig.1. This contour is also known as theKeldysh contour [20, 21]. A more detailed inspection ofEq.(4) then shows that the expectation value can also bewritten as a contour-ordered product [20, 22, 23, 24, 25].The one-particle Green function is then defined as acountour-ordered product of a creation and an annihi-lation operator G (1 ,
2) = − i h T C [ ˆ ψ H (1) ˆ ψ † H (2)] i , (5)where T C denotes the time-ordering operator on the con-tour and where we used the compact notation 1 = ( x , t )and 2 = ( x , t ). If we consider the Green function attime t = t − iβ and use the cyclic property of the trace,we find that G ( x t − iβ,
2) = − G ( x t ,
2) [24]. Hence,
FIG. 1: Keldysh contour. The depicted contour allows forthe calculation of observables for times t ≤ t ≤ T . Theinitial Green function is calculated on the imaginary track[ t , t − iβ ]. As we propagate the KB equations in time, forreal times t > t , the turning point of the time-contour at t = T moves to the right along the real time axis. the Green function defined in Eq. (5) obeys the boundaryconditions G ( x t ,
2) = − G ( x t − iβ, , (6) G (1 , x t ) = − G (1 , x t − iβ ) . (7)The Green function satisfies the equation of motion[ i∂ t − h (1)] G (1 ,
2) = δ (1 ,
2) + Z C d , G (3 , , (8)as well as a corresponding adjoint equation [9, 24]. InEq.(8) the time-integration is carried out along the con-tour C . The self-energy Σ incorporates the effects of ex-change and correlation in many-particle systems and is afunctional of the Green function that can be defined dia-grammatically [9, 19]. The Green function can be writtenas G (1 ,
2) = θ ( t, t ′ ) G > (1 ,
2) + θ ( t ′ , t ) G < (1 , , (9)where θ is a step function generalized to arguments on thecontour i.e. with θ ( t, t ′ ) = 1 if t is later on the contourthan t ′ and zero otherwise [9]. The greater and lessercomponents G > and G < respectively, have the explicitform G > (1 ,
2) = − i h ˆ ψ H (1) ˆ ψ † H (2) i , (10) G < (1 ,
2) = i h ˆ ψ † H (2) ˆ ψ H (1) i . (11)When one of the arguments is on the vertical track of thecontour, we adopt the notation [22] G ⌉ (1 , x , − iτ ) = G < (1 , x , t − iτ ) , (12) G ⌈ ( x , − iτ ,
2) = G > ( x , t − iτ , . (13)Finally, for the case when both time arguments are onthe imaginary track of the contour, we have the so-calledMatsubara Green function iG M [19] iG M ( x τ , x τ ) = G ( x t − iτ , x t − iτ ) , (14)which is a well-known object from the equilibrium theory.The factor i in the definition of Eq.(14) is a conventionwhich ensures that G M is a real function. The self-energyΣ has a similar general structure as the Green functionΣ(1 ,
2) = Σ HF (1 ,
2) + θ ( t, t ′ )Σ > (1 ,
2) + θ ( t ′ , t )Σ < (1 , . (15)The main difference with Eq.(9) is the appearance of theterm Σ HF which is proportional to a contour delta func-tion δ ( t , t ) in the time coordinates [9]. This term hasthe explicit formΣ HF [ G ](1 ,
2) = δ ( t , t )Σ HF ( x , x , t ) , (16)where Σ HF ( x , x , t ) = iG < ( x t, x t ) v ( x , x ) − iδ ( x − x ) Z d x v ( x , x ) G < ( x t, x t ) . (17)The structure of this self-energy is that of the Hartree-Fock (HF) approximation. However, in general we willevaluate this expression for Green functions G obtainedbeyond HF level (see section III). Using the form ofthe self-energy of Eq.(15) the contour integrations canbe readily carried out [9, 20] and we find separate equa-tions for the different Green functions G ≶ , G ⌉⌈ and G M .To display their temporal structure more clearly we sur-press the spatial indices of the Green functions and self-energies. Alternatively, these quantities may be regardedas matrices [7]. On the imaginary track of the contourwe obtain[ − ∂ τ − h ] G M ( τ ) = δ ( τ ) + Z β d ¯ τ Σ M ( τ − ¯ τ ) G M (¯ τ ) , (18)where the Green function and the self-energy are func-tions of the time-differences only, i.e. iG M ( τ − τ ) = G ( − iτ , − iτ ) and i Σ M ( τ − τ ) = Σ( − iτ , − iτ ), sincethe Hamiltonian is time-independent (and equal to ˆ H )on the imaginary track. Equation (18), which determinesthe Green function of the equilibrium system, has beentreated in detail in references [15, 16]. For the otherGreen functions we obtain i∂ t G ≶ ( t, t ′ ) = h HF ( t ) G ≶ ( t, t ′ ) + I ≶ ( t, t ′ ) , (19) − i∂ t ′ G ≶ ( t, t ′ ) = G ≶ ( t, t ′ ) h HF ( t ′ ) + I ≶ ( t, t ′ ) , (20) i∂ t G ⌉ ( t, − iτ ) = h HF ( t ) G ⌉ ( t, − iτ ) + I ⌉ ( t, − iτ ) , (21) − i∂ t G ⌈ ( − iτ, t ) = G ⌈ ( − iτ, t ) h HF ( t ) + I ⌈ ( − iτ, t ) , (22)where h HF ( t ) = h ( t ) + Σ HF ( t ) and Σ HF ( t ) is given byEq.(17). The retarded and advanced functions for G andΣ are defined according to F R/A ( t, t ′ ) = ± θ ( ± t ∓ t ′ )[ F > ( t, t ′ ) − F < ( t, t ′ )] , (23) with F replaced by G and Σ respectively. The so-calledcollision terms I ≶ and I ⌉⌈ have the form I ≶ ( t, t ′ ) = Z t d ¯ t Σ R ( t, ¯ t ) G ≶ (¯ t, t ′ ) + Z t ′ d ¯ t Σ ≶ ( t, ¯ t ) G A (¯ t, t ′ )+ 1 i Z β d ¯ τ Σ ⌉ ( t, − i ¯ τ ) G ⌈ ( − i ¯ τ , t ′ ) , (24) I ≶ ( t, t ′ ) = Z t d ¯ tG R ( t, ¯ t )Σ ≶ (¯ t, t ′ ) + Z t ′ d ¯ tG ≶ ( t, ¯ t )Σ A (¯ t, t ′ )+ 1 i Z β d ¯ τ G ⌉ ( t, − i ¯ τ )Σ ⌈ ( − i ¯ τ , t ′ ) , (25) I ⌉ ( t, − iτ ) = Z t d ¯ t Σ R ( t, ¯ t ) G ⌉ (¯ t, − iτ )+ Z β d ¯ τ Σ ⌉ ( t, − i ¯ τ ) G M (¯ τ − τ ) , (26) I ⌈ ( − iτ, t ) = Z t d ¯ t G ⌈ ( − iτ, ¯ t )Σ A (¯ t, t )+ Z β d ¯ τ G M ( τ − ¯ τ )Σ ⌈ ( − i ¯ τ , t ) . (27)These equations are readily derived using the conversiontable of Ref. [20]. From the symmetry relations G ≶ ( t, t ′ ) † = − G ≶ ( t ′ , t ) , (28)Σ ≶ ( t, t ′ ) † = − Σ ≶ ( t ′ , t ) , (29)it follows that we only need to calculate G > ( t, t ′ ) andΣ > ( t, t ′ ) for t > t ′ and G < ( t, t ′ ) and Σ < ( t, t ′ ) for t ≤ t ′ .These equations imply that I ≶ , ( t, t ′ ) = − I ≶ , ( t ′ , t ) † . Wefurther have G ⌈ ( − iτ, t ) = G ⌉ ( t, − i ( β − τ )) † , (30)Σ ⌈ ( − iτ, t ) = Σ ⌉ ( t, − i ( β − τ )) † . (31)The symmetry relations (28) and (30) for the Green func-tion follow directly from its definition, whereas the sym-metry relations (29) and (31) for the self-energy followfrom Eqs.(3.19) and (3.20) of Ref. [9]. Another conse-quence of equations (30) and (31) is that I ⌈ ( − iτ, t ) = (cid:2) I ⌉ ( t, − i ( β − τ )) (cid:3) † , which means that in practice it issufficient to calculate only I > , I < and I ⌉ . Eqs.(19) to(22) are known as the Kadanoff-Baym equations [8, 9].Once the Matsubara Green function G M ( τ ) is obtainedfrom Eq.(18), the Green functions G x ( x = ≶ , ⌉⌈ ) can becalculated by time propagation. Their initial conditionsare G > (0 ,
0) = iG M (0 + ) , (32) G < (0 ,
0) = iG M (0 − ) , (33) G ⌉ (0 , − iτ ) = iG M ( − τ ) , (34) G ⌈ ( − iτ,
0) = iG M ( τ ) . (35)The KB equations, together with the initial conditions,completely determine the Green functions for all timesonce a choice for the self-energy has been made. The formof the self-energy will be the topic of the next section. III. SELF-ENERGY APPROXIMATIONS
In the applications of the KB equations it is possible toguarantee that the macroscopic conservation laws, suchas those of particle, momentum and energy conservation,are obeyed. Baym [6] has shown that this is the casewhenever the self-energy is obtained from a functionalΦ[ G ], such that Σ(1 ,
2) = δ Φ δG (2 , . (36)Such approximations to the self-energy are called con-serving or Φ-derivable approximations. Well-known con-serving approximations are the Hartree-Fock, the secondBorn [8], the GW [26], and the T -matrix [8] approxi-mation. In our work we implemented the first three ofthese. The second Born approximation –
This approximationfor the self-energy consists of the two diagrams to secondorder in the two-particle interaction [8, 27]Σ(1 ,
2) = Σ HF (1 ,
2) + Σ (2) (1 , , (37)where Σ HF is the HF part of the self energy of Eq.(16)and Σ (2) = Σ (2 a ) + Σ (2 b ) is the sum of the two termsΣ (2 a ) (1 ,
2) = − i G (1 , Z d d v (1 , × G (3 , G (4 , v (4 , , (38)Σ (2 b ) (1 ,
2) = i Z d d G (1 , v (1 , G (3 , × G (4 , v (3 , , (39)where v (1 ,
2) = v ( x , x ) δ ( t , t ). These terms are usu-ally referred to as the second-order direct and exchangeterms. This approximation to the self-energy has beendiscussed in detail for the equilibrium case in Ref. [15].For the nonequilibrium case we need to calculate the vari-ous components Σ x ( x = ≶ , ⌉⌈ ). These are explicitly givenby Σ (2 a ) , ≶ (1 ,
2) = − i G ≶ (1 , Z d d v (1 , × G ≶ (3 , G ≷ (4 , v (4 , , (40)Σ (2 a ) , ⌉⌈ (1 ,
2) = − i Z d d G ⌉⌈ (1 , v (1 , × G ⌉⌈ (3 , G ⌈⌉ (4 , v (4 , , (41)for the direct diagram, andΣ (2 b ) , ≶ (1 ,
2) = i Z d d G ≶ (1 , v (1 , G ≷ (3 , × G ≶ (4 , v (3 , , (42)Σ (2 b ) , ⌉⌈ (1 ,
2) = i Z d d G ⌉⌈ (1 , v (1 , G ⌈⌉ (3 , × G ⌉⌈ (4 , v (3 , , (43) for the second-order exchange diagram. These expres-sions follow immediately from Eqs.(38) and (39) withhelp of the conversion table of Ref. [20]. The GW approximation – In the GW approximation theexchange-correlation part of the self-energy is given asa product of the Green function G with a dynamicallyscreened interaction W [26] . The screened interaction W satisfies the equation W (1 ,
2) = v (1 ,
2) + Z d d v (1 , P (3 , W (4 , . (44)Here, v is the bare Coulomb interaction, and P (1 ,
2) = − iG (1 , G (2 , , (45)is the irreducible polarization [26]. However, since thefirst term in Eq.(44) is singular in time (proportional toa delta function) it is convenient, for numerical purposes,to define its time-nonlocal part ˜ W = W − v [16]. FromEq.(44) it follows that˜ W (1 ,
2) = Z d d v (1 , P (3 , v (4 , Z d d v (1 , P (3 ,
4) ˜ W (4 , . (46)In terms of ˜ W , the self-energy has the form [26]Σ(1 ,
2) = Σ HF (1 ,
2) + iG (1 ,
2) ˜ W (1 , . (47)The part Σ c = iG ˜ W represents the correlation part ofthe self-energy and has the componentsΣ ≶ c (1 ,
2) = iG ≶ (1 ,
2) ˜ W ≷ (2 , , (48)Σ ⌉⌈ c (1 ,
2) = iG ⌉⌈ (1 ,
2) ˜ W ⌈⌉ (2 , . (49)From the fact that ˜ W (1 ,
2) has the same symmetriesas the contour-ordered density response function [26] χ (1 ,
2) = − i h T C [ˆ n H (1)ˆ n H (2)] i , where ˆ n is the densityoperator, it follows that˜ W ≶ (2 ,
1) = ˜ W ≷ (1 ,
2) = − [ ˜ W ≷ (2 , ∗ , (50) W ⌉ (1 ,
2) = W ⌈ (2 , . (51)In the following, we will again surpress the spatial co-ordinates in order to display the temporal structure ofthe equations more clearly. From the symmetry relations(50), (51), (48) and (49), and the fact that we only needΣ > ( t, t ′ ) for t > t ′ and Σ < ( t, t ′ ) for t ≤ t ′ , it follows thatwe only need to calculate ˜ W ⌉ ( t, − iτ ), and ˜ W < ( t, t ′ ) for t ≤ t ′ . The latter obey the equations:˜ W < ( t, t ′ ) = vP < ( t, t ′ ) v + vX < ( t, t ′ ) , (52)˜ W ⌉ ( t, − iτ ) = vP ⌉ ( t, − iτ ) v + vX ⌉ ( t, − iτ ) , (53)where P < ( t, t ′ ) = − iG < ( t, t ′ ) G > ( t ′ , t ) , (54) P ⌉ ( t, − iτ ) = − iG ⌉ ( t, − iτ ) G ⌈ ( − iτ, t ) , (55)and where the terms X < and X ⌉ are given by X < ( t, t ′ ) = Z t ′ d ¯ tP < ( t, ¯ t ) ˜ W A (¯ t, t ′ )+ Z t d ¯ tP R ( t, ¯ t ) ˜ W < (¯ t, t ′ )+ Z β d ¯ τ P ⌉ ( t, − i ¯ τ ) ˜ W ⌈ ( − i ¯ τ , t ′ ) , (56) X ⌉ ( t, − iτ ) = Z t ′ d ¯ tP R ( t, ¯ t ) ˜ W ⌉ (¯ t, − iτ )+ Z β d ¯ τ P ⌉ ( t, − i ¯ τ ) ˜ W M (¯ τ − τ ) , (57)with the retarded and advanced quantities defined as inEq.(23). The initial conditions for ˜ W < and ˜ W ⌉ are˜ W < (0 ,
0) = i ˜ W M (0 − ) , (58)˜ W ⌉ (0 , − iτ ) = i ˜ W M ( − τ ) , (59)where i ˜ W M ( τ − τ ′ ) = ˜ W ( t − iτ, t − iτ ′ ) is the Matsubarainteraction discussed in detail in Ref. [16]. IV. TIME-PROPAGATION OF THEKADANOFF-BAYM EQUATIONS
In the following, we will describe the time-propagationmethod which we employed to solve the KB equations.This method can be applied to general Hamiltonians con-taining one- and two-body interactions, and is further in-dependent of the explicit form of the self-energy.The time-propagation method is applied to the KBequations in matrix form. This matrix form is obtainedby expressing the Green function in terms of a set of ba-sis functions φ i ( x ), which we choose to be Hartree-Fockorbitals [7, 15, 16] G ( x t, x ′ t ′ ) = X ij G ij ( t, t ′ ) φ i ( x ) φ ∗ j ( x ′ ) . (60)When Eq.(60) is inserted in the expressions for the self-energy we obtain a basis set representation of the self-energy involving the matrices G ij ( t, t ′ ) and the two-electron integrals which are given as integrals of orbitalproducts with the two-body interaction v . All the quan-tities therefore become time-dependent matrices and allproducts are to be interpreted as matrix products. Wewill, however, surpress all matrix indices to display thetemporal structure of the equations more clearly. Ex-plicit expressions of the matrix form of the second Bornand GW self-energy are given in Refs.[7, 15, 16].We start by discussing the time-propagation of G > and G < . Due to the symmetry relations Eq.(28) and (29) weonly need to calculate G > ( t, t ′ ) for t > t ′ and G < ( t, t ′ )for t ≤ t ′ . From Eqs.(19) and (20) it then follows that G > must be time-stepped in the first time-argument FIG. 2: Time-stepping in the ( t, t ′ )-plane. G > ( t, t ′ ) is calcu-lated for t > t ′ and G < ( t, t ′ ) is calculated for t ≤ t ′ . and G < in the second one. We thus need to calcu-late G > ( T + ∆ , t ′ ) and G < ( t, T + ∆) for a small timestep ∆, from the knowledge of G ≷ ( t, t ′ ) for t, t ′ ≤ T .The symmetry relations (28) then immediately provideus with G > ( t ′ , T + ∆) and G < ( T + ∆ , t ) as well. Thetime-stepping procedure is illustrated in Fig.2 that dis-plays the ( t, t ′ )-plane, in which at a given time T all thequantities inside the square with sides equal to T , areknown. The time-step G < ( t, T ) → G < ( t, T + ∆) corre-sponds to a shift of the upper side of the time square with∆ i.e. a shift from the solid to the dotted line in Fig.2.Similarly the time-step G > ( T, t ′ ) → G > ( T + ∆ , t ′ ) corre-sponds to a shift of the righthand side of the time squarewith ∆. We further need to make a step G < ( T, T ) → G < ( T + ∆ , T + ∆) along the time diagonal t = t ′ . Thepropagation of G ⌈ ( − iτ, t ) and G ⌉ ( t, − iτ ) requires a time-step in the real time coordinate t for fixed imaginary timepoints τ .Note that the righthand sides of Eqs.(19) to (22) dependon the Green functions at the times T + ∆, which are notknown at time T . We therefore carry out the time-step T → T + ∆ twice. After taking the time step for the firsttime, we recalculate the righthand sides of Eqs.(19) to(22) and repeat the time-step T → T + ∆ using an aver-age of the old and new collision and HF terms. Since theterm h HF ( t ) in Eqs.(19) to (22) can attain large values, itis favorable to eliminate this term from the time-steppingequations. For each time-step T → T + ∆ we thereforeabsorb the term in a time-evolution operator of the form U ( t ) = e − i ¯ h HF ( T ) t , (61)where ¯ h HF ( T ) = h ( T + ∆ /
2) + Σ HF ( T ), where h is theone-body part of the Hamiltonian of Eq.(2). The one-body Hamiltonian h ( t ) is explicitly known as a functionof time and can be evaluated at half the time-step. Theterm Σ HF is only known at time T and will be recalcu-lated in the repeated time-step. In terms of the opera-tor U ( t ) of (61) we define new Green function matrices g x ( x = ≶ , ⌉⌈ ), as G ≶ ( t , t ) = U ( t ) g ≶ ( t , t ) U † ( t ) , (62) G ⌉ ( t , − iτ ) = U ( t ) g ⌉ ( t , − iτ ) , (63) G ⌈ ( − iτ , t ) = g ⌉ ( − iτ , t ) U † ( t ) . (64)We can now transform Eqs.(19) to (22) into equationsfor g x . For instance, g > satisfies the equation i∂ t g > ( t, t ′ ) = U † ( t )( h HF ( t ) − ¯ h HF ) G > ( t, t ′ ) U ( t ′ )+ U † ( t ) I > ( t, t ′ ) U ( t ′ ) . (65)Since ¯ h HF ≈ h HF ( t ) for times T ≤ t ≤ T + ∆, we canneglect for these times the first term on the right handside of Eq.(65). We then find G > ( T + ∆ , t ) = U ( T + ∆) g > ( T + ∆ , t ) U † ( t ) == U ( T + ∆) " g > ( T, t ) + Z T +∆ T dt∂ t g > ( t, t ) U † ( t ) ≈ U (∆) G > ( T, t ) − iU (∆) (Z T +∆ T dt e i ¯ h HF ( t − T ) ) I > ( T, t )= U (∆) G > ( T, t ) − V (∆) I > ( T, t ) , , (66)where V (∆) is defined as V (∆) = (¯ h HF ) − [1 − e − i ¯ h HF ∆ ] . (67)Similarly for G < , which is propagated using Eq.(20), wefind the equation G < ( t , T + ∆) = G < ( t , T ) U † (∆) − I < ( t , T ) V † (∆) . (68)For time-stepping along the time-diagonal we use i∂ t G < ( t, t ) = [ h HF ( t ) , G < ( t, t )]+ I < ( t, t ) − I < ( t, t ) , (69)which follows directly from a combination of the equa-tions for G < of Eqs.(19) and (20). The correspondingequation for g < ( t, t ) on the time diagonal then becomes i∂ t g < ( t, t ) = U † ( t )[ h HF ( t ) − ¯ h HF , G < ( t, t )] U ( t )+ U † ( t )( I < ( t, t ) − I < ( t, t )) U ( t ) . (70)From this equation we then obtain G < ( T + ∆ , T + ∆) == U ( T + ∆) g < ( T + ∆ , T + ∆) U † ( T + ∆)= U (∆) G < ( T, T ) U † (∆) − iU (∆) "Z ∆0 dt U † ( t ) I U ( t ) U † (∆) , (71) where we defined I = I < ( T, T ) − I < ( T, T ). By usingthe operator expansion e A Be − A = B + [ A, B ] + 12 [ A, [ A, B ]]+ 13 [ A,
12 [ A, [ A, B ]]] + . . . , (72)it follows that − i Z ∆0 dt U † ( t ) I U ( t ) = ∞ X n =0 C ( n ) , (73)where C ( n ) = i ∆ n + 1 [¯ h HF , C ( n − ] , (74)and C (0) = − i ∆ I . If we insert Eq.(73) into Eq.(71) wefinally obtain G < ( T +∆ , T +∆) = U (∆) " G < ( T, T ) + ∞ X n =0 C ( n ) U † (∆)(75)We found that keeping terms for n ≤ g ⌉ we have the equation i∂ t g ⌉ ( t, − iτ ) = U ( t ) † ( h HF ( t ) − ¯ h HF ) G ⌉ ( t, − iτ )+ U ( t ) † I ⌉ ( t, − iτ ) . (76)This yields, similarly as in Eq.(66) and (68) G ⌉ ( T + ∆ , − iτ ) = U (∆) G ⌉ ( T, − iτ ) − V (∆) I ⌉ ( T, − iτ ) . (77)Finally, for G ⌈ we have G ⌈ ( − iτ , T + ∆) = G ⌈ ( − iτ , T ) U (∆) † − I ⌈ ( − iτ , T ) V (∆) † . (78)The Eqs. (66), (68), (71), (77) and (78) form the basisof the time-stepping algorithm. At each time-step, itrequires the construction of the step operators U (∆)and V (∆) and therefore the diagonalization of ¯ h HF forevery time-step. As mentioned before, the righthandsides of Eqs.(19) to (22) depend on the Green functionsat the times T + ∆ which are not known at time T . Wetherefore carry out the time-step T → T + ∆ twice. Theprocedure is as follows:(1) The collision integrals and ¯ h HF at time T arecalculated from the Green functions for times t, t ′ ≤ T .(2) A step in the Green function G ( T ) → G ( T + ∆) istaken according to Eqs.(66), (68), (71), (77) and (78).(3) New collision integrals I > ( T + ∆ , t ) , I > ( t, T +∆) , I ⌉ ( T + ∆ , − iτ ) and I ⌈ ( − iτ, T + ∆) are calculated byinserting the new Green functions for times t, t ′ ≤ T + ∆into Eqs.(24) to (27).(4) The values of the collision integrals andthe Hartree-Fock self-energy are approximatedby ¯ I = ( I ( T ) + I ( T + ∆)) / HF =(Σ HF ( T ) + Σ HF ( T + ∆)) / I ( T ) and I ( T + ∆)are the collision terms calculated under points (1) and(3).(5) The Green function is then propagated from G ( T ) → G ( T + ∆) using the average values ¯ I and¯ h HF = h ( T + ∆ /
2) + ¯Σ HF in Eqs.(66), (68), (71), (77)and (78).This concludes the general time-stepping procedurefor the Green functions.We finally consider the calculation of ˜ W < and ˜ W ⌉ fromEqs. (52) and (53). As a consequence of the symmetryrelation (50), we only need to calculate ˜ W < ( t, t ′ ) for t < t ′ . In a time step from T to T + ∆ we needto calculate ˜ W < ( t, T + ∆) for t ≤ T + ∆ from theknown values of ˜ W < ( t, T ) for t ≤ T . The first termon the righthand side of Eq.(52) can be calculateddirectly from G < ( t, T + ∆) and G > ( T + ∆ , t ). However,the last term X < ( t, T + ∆) of Eq.(52) depends onthe, still undetermined, values ˜ W < ( t, T + ∆). Wetherefore employ an iterative scheme. As a first guessfor ˜ W < ( t, T + ∆) we take ˜ W < ( t, T + ∆) = ˜ W < ( t, T )for t ≤ T and ˜ W < ( T + ∆ , T + ∆) = ˜ W < ( T, T ). Wetherefore use the values of ˜ W < on the known sides ofthe time square at time T (solid lines in Fig.2) as initialguesses for the stepped sides (dotted lines in Fig.2) at T + ∆. As an initial guess for the value of ˜ W < at thenew diagonal point ( T + ∆ , T + ∆), we take the value atthe previous diagonal point ( T, T ). We then calculatethe quantity X < ( t, T + ∆) for t ≤ T + ∆ and obtain anew value for ˜ W < ( t, T + ∆) from Eq.(52). This value is then inserted back into the righthand side of Eq.(52)and the process is repeated until convergence is reached.Similarly we initialize ˜ W ⌉ ( T + ∆ , − iτ ) = ˜ W ⌉ ( T, − iτ )and solve Eq.(53) in the same manner as for ˜ W < .This concludes our derivation of the time-stepping al-gorithm of the KB equations. The propagation methoddescribed here has been used in two recent Letters [7, 11]where also values for the numerical parameters aregiven. It is clear that the choice of these parametersdepends strongly on the type of system considered, andon the strength of the applied external fields. V. SUMMARY AND CONCLUSIONS
We presented a detailed account of the KB equationsand discussed in detail their structure, initial conditionsand symmetries. We developed an algorithm for thetime-propagation of the KB equations in which the sym-metry relations for the Green functions were used toreduce the set equations that needed to be solved. Intwo recent Letters [7, 11] we applied the method to thecase of atoms and molecules in external time-dependentfields and to the case of transient transport dynamics ofdouble quantum dots. We therefore conclude that time-propagation of the KB equations can be used as a practi-cal method to calculate the nonequilibrium properties ofa wide variety of many-body quantum systems, rangingfrom atoms and molecules to quantum dots and quan-tum wells. Moreover, the present work can be readilyextended to other Green function formalisms, such as theNambu formalism [28, 29] for superconducting systems.Also future extension to bosonic systems is straightfor-ward. Work along these lines is in progress. [1] G. Cuniberti, G. Fagas, and K. Richter, eds.,
Introduc-ing Molecular Electronics , vol. 680 (Springer, New York,2005).[2] T. H. T. Kato, T. Aoki, T. Ando, A. Fukuda, and S.-S.Seomun, Phys. Rev. Lett. , 226804 (2003).[3] J. M. Elzerman, R. Hanson, L. H. W. van Beveren,B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwen-hoven, Nature , 431 (2003).[4] J. Shah, ed., Ultrafast Spectroscopy of Semiconductorsand Semiconductor Nanostructures (Springer, Berlin,1999).[5] M. Merano, S. Sonderegger, A. Crottini, S. Collin,P. Renucci, E. Pelucchi, A. Malko, M. Baier, E. Kapon,B. Deveaud, et al., Nature , 479 (2007).[6] G. Baym, Phys. Rev. , 1391 (1962).[7] N. E. Dahlen and R. van Leeuwen, Phys. Rev. Lett. ,153004 (2007).[8] L. P. Kadanoff and G. Baym, Quantum Statistical Me-chanics (W. A. Benjamin, Inc., New York, 1962).[9] P. Danielewicz, Ann. Phys. (N. Y.) , 239 (1984). [10] N.-H. Kwong and M. Bonitz, Phys. Rev. Lett. , 1768(2000).[11] P. My¨ohanen, A. Stan, G. Stefanucci, and R. vanLeeuwen, Europhys. Lett. , 67001 (2008).[12] W. Sch¨afer, J.Opt.Soc.Am. B13 , 1291 (1996).[13] M. Bonitz, D. Kremp, D. C. Scott, R. Binder, W. Kraeft,and H. S. K¨ohler, J.Phys.Cond.Matter , 6057 (1996).[14] R. Binder, H. S. K¨ohler, M. Bonitz, and N. Kwong,Phys.Rev. B55 , 5110 (1997).[15] N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. ,164102 (2005).[16] A. Stan, N. E. Dahlen, and R. van Leeuwen, J. Chem.Phys. (accepted) (2009).[17] H. S. K¨ohler, N. H. Kwong, and H. A. Yousif, Comp.Phys. Comm. , 123 (1999).[18] M. Bonitz and D. Semkat, in
Introduction to Com-putational Methods in Many-Body Physics , edited byM. Bonitz and D. Semkat (Rinton Press, Princeton,2006), p. 171.[19] A. L. Fetter and J. D. Walecka,
Quantum Theory of
Many-Particle Systems (McGraw-Hill, New York, 1971).[20] R. van Leeuwen, N. E. Dahlen, G. Stefanucci, C.-O. Alm-bladh, and U. von Barth, in
Time-dependent DensityFunctional Theory , edited by M. A. L. Merques, C. A.Ullrich, F. Nogueira, A. Rubio, K. Burke, and E. K. U.Gross (Springer, Berlin Heidelberg, 2006), p. 33.[21] L. V. Keldysh, Zh. Eksp. Teor. Fiz. , 1515 (1964), [ Sov.Phys. JETP , , 1018 (1965)].[22] G. Stefanucci and C.-O. Almbladh, Phys. Rev. B ,195318 (2004).[23] M. Wagner, Phys. Rev. B , 6104 (1996). [24] N. E. Dahlen, A. Stan, and R. van Leeuwen, J. Phys.Conf. Ser. , 324 (2006).[25] N. E. Dahlen, R. van Leeuwen, and A. Stan, J. Phys.Conf. Ser. , 340 (2006).[26] L. Hedin, Phys. Rev. , A796 (1965).[27] G. Baym and L. P. Kadanoff, Phys. Rev. , 287 (1961).[28] J. R. Schrieffer, Theory of Superconductivity (Addison-Wesley, 1988).[29] Y. Nambu, Phys.Rev.117