The Kadanoff-Baym approach to double excitations in finite systems

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. The Kadanoff-Baym approach to double excitations in finite systems Säkkinen, Niko; Manninen, Matti; van Leeuwen, Robert Säkkinen, N., Manninen, M., & van Leeuwen, R. (2012). The Kadanoff-Baym approach to double excitations in finite systems. New Journal of Physics, 14(13032). https://doi.org/10.1088/1367-2630/14/1/013032 2012 This content has been downloaded from IOPscience. Please scroll down to see the full text. Download details: IP Address: 130.234.75.141 This content was downloaded on 14/01/2016 at 07:38 Please note that terms and conditions apply. The Kadanoff–Baym approach to double excitations in finite systems View the table of contents for this issue, or go to the journal homepage for more 2012 New J. Phys. 14 013032 (http://iopscience.iop.org/1367-2630/14/1/013032) Home Search Collections Journals About Contact us My IOPscience The open–access journal for physics New Journal of Physics The Kadanoff–Baym approach to double excitations in finite systems N S¨ akkinen1, M Manninen and R van Leeuwen Department of Physics, Nanoscience Center, University of Jyv¨ askyl¨ a, FIN 40014 Jyv¨ askyl¨ a, Finland E-mail: [email protected] New Journal of Physics 14 (2012) 013032 (21pp) Received 18 August 2011 Published 18 January 2012 Online at http://www.njp.org/ doi:10.1088/1367-2630/14/1/013032 Abstract. We benchmark many-body perturbation theory by studying neutral, as well as non-neutral, excitations of finite lattice systems. The neutral excitation spectra are obtained by time-propagating the Kadanoff–Baym equations in the Hartree–Fock and the second Born approximations. Our method is equivalent to solving the Bethe–Salpeter equation with a high-level kernel while respecting self-consistency, which guarantees the fulfillment of a frequency sum rule. As a result, we find that a time-local method, such as Hartree–Fock, can give incomplete spectra, while already the second Born approximation, which is the simplest time-non-local approximation, reproduces well most of the additional excitations, which are characterized as double excitations. 1Author to whom any correspondence should be addressed. New Journal of Physics 14 (2012) 013032 1367-2630/12/013032+21$33.00 © IOP Publishing Ltd and Deutsche Physikalische Gesellschaft 2 Contents 1. Introduction 2 2. Response functions in non-equilibrium many-body theory 3 2.1. The non-equilibrium Kadanoff–Baym approach ................. 3 2.2. Generalized response function .......................... 6 2.3. Approximate Bethe–Salpeter kernels ....................... 7 3. Model and test systems 8 4. Results: addition, removal and excitation spectra 10 4.1. Numerical aspects ................................. 10 4.2. Addition and removal energies .......................... 10 4.3. Excitation energy spectra ............................. 13 5. Summary and conclusions 17 Acknowledgments 18 Appendix. The Thomas–Reiche–Kuhn/f-sum rule 18 References 20 1. Introduction We recently developed a non-equilibrium many-body method, based on time propagation of the Kadanoff–Baym equations, to study transient dynamics in quantum transport, and found that electronic correlations beyond mean field strongly affect the relevant non-equilibrium properties [1,2]. The method is approximate, because only some selected classes of perturbative terms are accounted for and therefore an estimate of the quality of the results would be desirable. However, as benchmark results are scarcely available for such open systems, we need to devise simpler experiments to test our method. Some finite, exactly solvable lattice systems, which typically model molecular devices and whose properties thus determine the transport characteristics, form an appropriate test platform. A many-body approximation which is unable to reproduce the properties of such an isolated system correctly has no prerequisites to fare well for a molecular device between conducting leads. The Kadanoff–Baym equations are equations of motion for the one-body Green’s function G, where many-body effects are incorporated into the so-called self-energy 6. The latter is a functional of the Green’s function, and the main task of many-body perturbation theory (MBPT) is to find good approximations for this quantity. We are, in particular, interested in the so-called conserving approximations for the self-energy, which are vital in quantum transport [3,4] as they guarantee the fulfillment of the conservation laws for the particle number, total momentum and energy [5,6]. Some of these approximations have already been tested in finite lattice systems [7,8]. These reports have highlighted strongly perturbed systems, in which serious issues, such as artificial damping of the dynamics, were observed and deemed correlation induced as they were not present in a mean-field approximation. The purpose of this paper is to extend these investigations to the linear response regime. Within this regime, the key quantity is the response function δG/δv which describes the change in the Green’s function due to a weak external field vand gives direct access to the neutral excitation spectra. This quantity satisfies an integral equation known as the New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 3 Bethe–Salpeter equation (BSE) [9] whose integral kernel is given by the functional derivative δ6/δG. Using this equation to calculate the excitation spectra has proven to be computationally challenging and, in practice, often requires a number of additional approximations, such as neglect of self-consistency, kernel diagrams and/or frequency dependence. Instead, we obtain the response function by time propagation of the Green’s function which does not require any of the above-mentioned approximations [10,11]. Moreover, such excitation spectra obtained automatically satisfy a frequency sum rule known as the f-sum rule, which acts as an important consistency relation for the response function. The spectral properties of neutral, as well as non-neutral, excitations depend strongly on the temporal structure of the self-energy which can be either time local or time non-local. Whereas a time-local scheme leads to merely renormalized one-particle excitations, a time-non-local one can reproduce a more complicated spectral structure. Such excitations involving double- or many-particle character are, in fact, dominant for some potential materials for realistic devices, such as organic polymers (polyenes) [12], making time-non-local approximations an interesting research topic. This represents a major challenge for the Bethe–Salpeter approach as time non-locality translates to an integral kernel which depends on three frequencies. Although some state-of-the-art calculations have recently been carried out in finite systems with a timedependent non-local kernel [13,14], these results have come at the expense of self-consistency, which is a requirement for the fulfillment of conservation laws and is therefore indispensable in a time-dependent transport process. Moreover, a non-self-consistent approximation can lead to a gross violation of the f-sum rule [10,15]. On the other hand, our calculations, performed within the time-local Hartree–Fock (HF) and time-non-local second Born (2B) approximations [11,16], do not sacrifice self-consistency, remain computationally tractable and yield excitation spectra satisfying the f-sum rule. The paper is organized as follows. We start with an introductory section 2in which we summarize the non-equilibrium many-body formalism, describe our means to approach the neutral excitation spectra and introduce, as well as analyze, our many-body approximations. In section 3, we introduce the test environment, which comprises two simple lattice systems. Our numerical results, in particular on addition, removal and excitation spectra, as well as the related discussion are contained in section 4. We summarize our work and present the conclusions in section 5. Finally, the appendix contains an extension of the proof of the f-sum rule to lattice systems. 2. Response functions in non-equilibrium many-body theory 2.1. The non-equilibrium Kadanoff–Baym approach The Green’s functions are reduced quantities that act as probes to physical observables while containing only a minimal amount of excess information. The non-equilibrium one-body Green’s function is defined as Gi j (z,z0)≡ −ihTγˆcH,i(z)ˆc† H,j(z0)i,(1) where h. . .idenotes an ensemble average with respect to the grand canonical density operator and ˆc(†) H,i(z)denotes Heisenberg picture operator which annihilates (creates) a particle from (to) a single-particle quantum state i. The contour-ordering operator Tγarranges the contour times zand z0along the extended Keldysh contour γ[17] as shown in figure 1; for further details, see [18–21]. New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 4 Figure 1. The extended Keldysh contour starting from t0and ending at t0−iβ consists of an equilibrium, imaginary track and two non-equilibrium, realtime branches. The zero-temperature limit β→ ∞, where βis the inverse temperature, is implied when finite systems are considered. The contour-ordered Green’s function, although a powerful formal device, does not represent directly any observable. It can, however, be reduced to a set of more intuitive and physical objects by exposing all possible time orderings. The one-body Green’s function, being a two-time object, can then be written as G(z,z0)=θ(z,z0)G>(z,z0)+θ(z0,z)G<(z,z0), (2) where the boldfaced italic symbols denote matrices in the basis spanned by the single-particle quantum states. The greater and lesser Keldysh components, which are defined as G> i j (z,z0)≡ −ihˆcH,i(z)ˆc† H,j(z0)i,(3a) G< i j (z,z0)≡ihˆc† H,j(z0)ˆcH,i(z)i,(3b) describe the motion of a particle and a hole, respectively, in the many-particle system. The core of the non-equilibrium theory lies in its equations of motion, the first of which can be written as (i∂z−h(z))G(z,z0)=δ(z,z0)+Zγ d¯z6[G](z,¯z)G(¯z,z0), (4) while the second, adjoint equation can be obtained by differentiation with respect to the second time argument. The symbol h(z)denotes the one-body part of the Hamiltonian, which can, in general, be time dependent. Using the Langreth rules [19,22], the equations of motion can be split into a set of equations for the different Keldysh components referred to as the Kadanoff–Baym equations [18,23]. The self-energy 6is an effective, non-local potential, which accounts for all the many-body effects, and allows a systematic, beyond order-by-order perturbation expansion. A Green’s function obtained from an approximate self-energy 6i j (z,z0)=δ8[G] δGji (z0,z),(5) where 8[G] is the Baym functional [6], obeys the conservation laws for the particle number, total momentum and energy provided it satisfies the equation of motion (4), or in other words is solved self-consistently [5,6]. The non-local structure of the one-body Green’s function ensures that all one-body observables which are accessible in the equilibrium/zero-temperature theory can now be calculated in non-equilibrium settings from the knowledge of the greater and lesser components. New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 5 In particular, all time-local observables are related to the one-body reduced density matrix (1RDM), which is given by γ(t)≡ −iG<(t,t), (6) where tis a real-time argument located on the horizontal track of the Keldysh contour. On the other hand, time non-locality allows access to information on particle number changing processes contained in the spectral function, which is defined as A(t,t0)≡iG>(t,t0)−G<(t,t0)(7) and whose pole structure in the frequency domain gives the particle addition and removal energies. Another virtue of time non-locality is the possibility of accessing some timelocal two-body observables, including the total energy which can be written in terms of the Galitski–Migdal (GM) functional as E(t)= − i 2tr[(i∂t+h)G<(t,t0)]t0=t,(8) where ∂tis proportional to the identity matrix, and trA≡PiAii is the regular algebraic trace. However, the main advantage of the non-equilibrium theory addressed here is the additional possibility of calculating some time-non-local two-body observables from the sole knowledge of the one-body Green’s function. We focus on the real-time, retarded response function, which is defined as the functional derivative χR i j,kl(t,t0)≡δγi j (t) δvlk (t0)v=0 (9) with respect to an external perturbation vi j (t), and has the general structure of a retarded function given by χR i j,kl(t,t0)=θ(t−t0)(χ> i j,kl(t,t0)−χ< i j,kl(t,t0)),(10) where the greater and lesser components are defined, respectively, in terms of the 1RDM operator ˆγH,i j (t)≡ ˆc† H,j(t)ˆcH,i(t), as χ> i j,kl(t,t0)≡ −ih ˆγH,i j (t)ˆγH,kl (t0)i,(11a) χ< i j,kl(z,z0)≡ −ih ˆγH,kl (t0)ˆγH,i j (t)i.(11b) The retarded response function therefore describes the motion of a particle–hole pair, or in other words an excitation, in the many-particle system. This quantity is connected to the one-body Green’s function through the linear response relation, which can be written as δγi j (t)=X kl Z∞ t0 d¯ tχR i j,lk(t,¯ t)vkl (¯ t), (12) where δγi j (t)is the first-order variation of the 1RDM. The lowest order response of the system to an external perturbation is therefore uniquely determined by the real-time, retarded response function whose pole structure in the frequency domain gives the neutral excitation energies. In the non-equilibrium approach, both the 1RDM and perturbation are known and therefore the response function is readily obtained by a proper choice of perturbation and subsequent inversion of equation (12) [10,11]. New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 6 L iz lz jzkz = iz lz jzkz +KL iz lz jzkz Figure 2. The Bethe–Salpeter equation for the generalized response function whose particular time ordering Li j,kl(zz+,z0), known as the particle–hole propagator, yields the retarded response function. The dashed lines represent possible connections, while double lines denote dressed Green’s functions. 2.2. Generalized response function The quality of the one-body Green’s function and, because of equation (12), also of the retarded response function is clearly determined by the self-energy approximation. Here we write down an explicit relation between the self-energy and the retarded response function. We start by introducing the generalized, contour-ordered response function, which is defined as Li j,kl(zz0,z00)≡δGi j (z,z0) δvlk (z00)v=0 (13) and in terms of which the greater and lesser components of the real-time response function, defined in equation (11), are given by χ≷ i j,kl(t,t0)= −iLi j,kl(zz+,z0)z=t±,z0=t0 ∓ ,(14) where z+≡z+ηwith the limit η→0+implied after contour ordering. The contour times z=t± are located at the lower/upper horizontal branch of the Keldysh contour shown in figure 1, and correspond to a real time t. The generalized response function can be straightforwardly shown [9] to obey the so-called Bethe–Salpeter equation, which is given by Li j,kl(zz0,z00)=Gil (z,z00)Gkj (z00,z0) +X pqrs Z γ d(z1z2z3z4)Gip(z,z1)Gq j (z2,z0)Kpq,sr (z1z2,z4z3)Lrs,kl(z3z4,z00), (15) where we defined the four-point Bethe–Salpeter kernel as the functional derivative Ki j,kl(zz0,¯z¯z0)≡δ6i j (z,z0) δGlk(¯z0,¯z).(16) Figure 2contains a diagrammatic representation of this equation, and illustrates the role of the kernel as a source of an infinite perturbation series. We can conclude that any retarded response function obtained by inversion of equation (12), and relating to a one-body Green’s function for an approximate self-energy 6, corresponds to a solution of the Bethe–Salpeter equation with the kernel δ6/δG. The correspondence further implies that a non-equilibrium approach already with a relatively simple self-energy accounts for quite a sophisticated approximation for the response functions (figure 3). New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 7 ΣHF = Hartree (H) + Fock (F) KHF = H + F (a) Σ2B= + + bubble (B) + exchange (X) K2B= + + + + B1−B3 + + + X1−X3 ( b ) Figure 3. The Hartree–Fock and second Born self-energies (6HF/2B) and Bethe–Salpeter kernels (KHF/2B =δ6HF/2B/δGHF/2B). Each self-energy diagram shown relates to at least one kernel diagram, for example, bubble (B) gives rise to diagrams B1−B3. Wiggly lines denote the Coulomb interaction, and double lines denote dressed Green’s functions. 2.3. Approximate Bethe–Salpeter kernels The time propagation of a Green’s function with an approximate self-energy 6automatically leads to a retarded response function related to an approximate kernel δ6/δG, which, for example, guarantees the obedience of the f-sum rule, as proven in the appendix. Such a consistency ensures that each self-energy diagram containing nGreen’s function lines gives rise to ninequivalent kernel diagrams which contribute to the response function. Our two approximate kernels are given by functional derivatives of the HF and 2B self-energies, and represent a time-local mean-field and a non-local correlated approximation, respectively. We can view our generalized response functions diagrammatically by inserting the approximate kernels into the Bethe–Salpeter equation shown in figure 2. Guided by this recipe, we easily deduce that the HF kernel shown in figure 4(a) yields a response function given by a simple bubbles-and-ladders series, whereas the 2B kernel, which is given in figure 4(b), results in a considerably more complicated response function. Note that, our approximations, whether related to a self-energy or a kernel, are known by the name of the self-energy without any extensions. For example, our HF approximation for the kernel shown in figure 4(a) is often explicitly referred to as the time-dependent HF approximation. We can also assign a meaning to each diagram in these perturbation expansions by considering a particular time ordering, which, in our case, is the ordering given by the particle–hole propagator Li j,kl(zz+,z0). We interpret our kernel diagrams as follows. The firstorder kernel diagrams denoted by H and F represent either a direct process (F) in which a particle–hole pair enters, interacts and exits the event or an exchange process (H) in which a particle–hole pair enters, annihilates and creates another pair, which then exits the event. The second-order kernel diagrams B1−3and X1−3can be deciphered similarly. Diagram B1describes a partially screened interaction between a particle and a hole, while diagrams B2–B3reflect the New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 14 Table 2. The exact ground state commutator c1≡c1 ED;11, and relative errors c1 HF/2B;11/c1−1, where c1 HF/2B;11 is the ground state commutator in the HF/2B approximation. U c1HF (%) 2B (%) (a) Half-filled, six-site, PPP chain 1.0 1.755 1.1 0.4 2.0 1.733 3.9 2.1 (b) One-sixth-filled, six-site, Hubbard ring 1.0 1.324 0.7 −0.2 2.0 1.307 2.1 1.2 where c1,(t/ω) i j ≡Pαβ c1,(t/ω) iα, jβ. On the left-hand side of this equation, we have a commutator (see equation (A.14)) defined as c1 i j ≡hji γi j +γji hi j −δi j X k (hjkγkj +γjk hkj ),(24) where hi j and γi j ≡Pαβ γiα, jβare the ground state, one-body Hamiltonian and reduced density matrix, respectively. In contrast, on the right-hand side, we have the first time derivative and frequency moment of the density response function, which are given, respectively, by c1,t i j ≡ −∂tχi j (t)|0+,(25a) c1,ω i j ≡ − 1 πZ∞ −∞ dω ω =m(χi j (ω)).(25b) We have also checked numerically, by fixing a method (ED, HF, 2B) and comparing the commutator (24) against the time derivative (25a), that the f-sum rule is satisfied up to a numerical accuracy better than 0.01% for all of our methods (ED, HF, 2B). Moreover, since our response functions are obtained by finite-length time propagations, we have reserved the comparison between the commutator (24) and frequency moment (25b) as a test of the completeness of our excitation spectra. As the frequency moment should be equal to the commutator, we can use the latter to estimate how well an approximate method is able to reproduce the true spectra. In table 2, we compare exact (ED) and approximate (HF, 2B) ground state commutators (24). Such a comparison shows that the 2B captures, in all cases, more precisely the value of the exact commutator, and is therefore expected to reproduce better-quality spectra than the HF. (iii) Lastly, we need some means of characterizing excitations which appear in the exact excitation spectra. As the expectation value of an occupation number operator contains information about single-particle state occupations in a compact statistical format, it is a suitable tool for the task. Our occupation number operator is given by ˆnHF i≡X α ˆ d† iαˆ diα,(26) where ˆ d(†) iαannihilates (creates) an electron from (to) the single-particle eigenstate {i, α} of the HF Hamiltonian. We can obtain information about single-particle transitions by New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 15 0 3 6 0 1 2 3 4 E1 E2 E3 E4 E5 ω -Imχ11(ω)/π U = 2.0 0 3 6 E1 E2 E3 E4 E5 0 1 2 3 4 U = 1.0 ED HF 2B Figure 8. The imaginary part of the density response function χ11(ω) of a sixsite/electron (half-filled) chain with long-range (PPP) interactions, for which the f-sum rule is satisfied up to 99.9–100.2% when c1,ω 11 is used. Exact excitation energies with a non-zero oscillator strength (ED-D) are given by the dashed, vertical lines (v1=0.01,T=150, 1 =0.02–0.025). comparing the expectation values of this operator for the ground and excited states of the many-particle system. However, since excitation spectra can comprise degenerate excited states which have varying spectral intensities, we need a construct that can distinguish the character of an excitation in such a situation. For this purpose, we define the weighted occupation numbers NHF i j,k(E)≡Tr(ρi j (E)ˆnHF k),(27) where Tr is the trace over a complete set of many-particle states, and ˆρi j (E)≡Pk,Ek=EPαβ fk,iαfk,jβ|8N kih8N k| Pk,Ek=EPαβ fk,iαfk,jβ ,(28) a properly normalized density operator which projects any state to a degenerate subspace, with eigenvalue E, of the many-body eigenstates. The oscillator strength fk,iαis defined below equation (21). In practice, we calculate the differences 1NHF i j,k(E)≡NHF i j,k(E)−NHF i j,k(E0), (29) where E0is the true ground state energy. This quantity tells us what kind of single-particle transitions are involved in the many-particle transitions from the ground state to the excited states with energy E. For example, consider a spin-compensated two-particle and two-level system with nearly one-determinental ground state. Then weighted occupation numbers for the lowest single-particle state whose values are close to one or two describe one- or twoparticle excitations, respectively. The neutral excitation spectrum of the half-filled chain is shown in figure 8for two (U= 1.0, 2.0) interaction strengths. The low interaction spectra (see the upper panel) are very similar; that is, both many-body approximations reproduce the intensities and frequency structure of the New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 16 -2 -1 0 1 2123456 ΔNk HF(E) E1 123456 E2 123456 E3 123456 E4 -2 -1 0 1 2 123456 E5 (a) -2 -1 0 1 2123456 ΔNk HF(E) E1 123456 E2 123456 E3 123456 E4 -2 -1 0 1 2 123456 E5 ( b ) Figure 9. The relative weighted occupation numbers, 1Nk(E)≡1NHF 11,k(E)for each Hartree–Fock state (k=1. . . 6) and excitation E1–E5 are shown for the half-filled chain (a) and one-sixth-filled ring (b) at U=1.0. The histograms whose filled (positive) or depleted (negative) parts sum up to 0.5–1.5 or 1.5–2.0 describe a state with a singly or doubly excited character, respectively. exact response function quite accurately, although the 2B approximation is slightly better. As the interaction strength is increased (see the lower panel) both many-body approximations still reproduce quantitatively correct results, although some low-energy structure, see excitations E1 and E2, has started to shift toward the lower end of the spectrum. This shift is a possible sign of a failure to describe the true ground state, as backed up by the fact that the squared overlap of the HF ground state with the true ground state is only 0.86. However, since the correction EGM–E0≈0.06 is small, such a shift cannot be attributed solely to a breakdown of the ground state. Overall, dominant features, such as the intensities and energies of low-energy excitations, as well as finer details, such as the diminishing intensity of the third dominant pole (E3) of the response function as a function of the interaction, are replicated more accurately by the correlated approximation, while the mean-field picture even fails to capture the subtleties. The fact that even the HF method reproduces all five, dominant excitation energies (E1–E5) correctly implies that we are dealing with one-particle excitations, as verified by the weighted occupation numbers of figure 9(a). We will shortly consider the question of how the manybody approximations cope with a more complicated excitation spectrum, in which additional excitations of many-particle nature arise. We give an example of such a system by introducing a less rigid, but energetically more simple, one-sixth-filled, six-site ring whose neutral excitation spectrum is given in figure 9for two (U=1.0, 2.0) interaction strengths. The results portray, irrespective of the interaction, adequate agreement between the exact and 2B excitation spectra, but, in contrast an utter failure of the HF approximation to account correctly for the structure of the response function. Two excitations, one at a lower energy (E2) and another at a higher energy (E4–E5), which are completely absent in HF, indicate the presence of two-particle excitations, which cannot be described with a time-local approximation. The weighted occupation numbers, shown in figure 9(b), confirm these speculations: only two excitations, E1 and E3, relate to singly excited states, whereas all others, in particular E2, E4 and E5, have a strong doubly excited state New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 17 0 2 4 0 1 2 3 4 5 6 E1 E2 E3 E4 E5 ω U = 2.0 -Imχ11(ω)/π 0 2 4 E1 E2 E3 E4 E5 0 1 2 3 4 5 6 U = 1.0 ED HF 2B Figure 10. The imaginary part of the density response function χ11(ω) of a sixsite, one-sixth-filled (two electron) ring with short-range (Hubbard) interactions, for which the f-sum rule is satisfied up to 99.8–100.1% when c1,ω 11 is used. Exact excitation energies with a non-zero oscillator strength (ED-D) are given by dashed, vertical lines (v1=0.01;T=150, 1 =0.02–0.025). character. Therefore, the highest excitation of HF attempts to mimic a non-existent excitation related dominantly to the lowest to highest state transition, see figure 6(b), whereas the 2B approximation successfully reproduces the two-particle excitations, the lowest one being quite accurate, whereas the higher ones are shifted toward lower energies. As a disadvantage of the correlated scheme, a possibly spurious excitation appears on the left-hand side of E3 when U=2.0. Such non-physical defects have also been seen in other recent studies [13,14], where they are associated with a response kernel which does not have a proper mathematical structure. 5. Summary and conclusions We have studied the performance of the time-dependent HF and 2B approximations by investigating the spectral properties of some finite, strongly correlated, lattice systems. Our purpose has been twofold: firstly, to test approximations used in transient quantum transport, and secondly, to gain insight into excitations in MBPT, both using the Kadanoff–Baym equations. We have calculated ground state energies, as well as spectral and density response functions with our approximate methods, and compared these against exact results, which were obtained by ED. The results show that both approximations perform well in a simple half-filled system. That is, they reproduce correctly the true excitation spectra, although the quasi-particle properties are not replicated as accurately, especially in the HF approximation. However, at a lower filling, we observe two-particle excitations, which by construction cannot be captured with the HF method, but are generated correctly in the 2B approximation. Such a difference is also seen in the spectral function in which additional quasi-particles are formed in this correlated approximation. Overall, the 2B approximation is consistently in better agreement with the exact New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 18 results, although signs of spurious excitations, also in more asymmetric and strongly interacting systems, call for further research. We conclude that a time-local approximation can be inadequate for highly correlated lattice systems at a low filling, which should have an impact on quantum transport, as particle number on device varies. In contrast, already the simplest time-non-local, self-consistent approximation performs more reliably. More complicated, but computationally realizable, many-body approximations, as well as investigations of the effects of self-consistency, remain subjects for a future work. Acknowledgments We acknowledge CSC—IT Center for Science Ltd for the allocation of computational resources and the Academy of Finland for financial support. Appendix. The Thomas–Reiche–Kuhn/f-sum rule The exact density(–density) response function satisfies a consistency relation known in atomic physics as the Thomas–Reiche–Kuhn (TRK) sum rule or in condensed matter physics as the f-sum rule. It is known [20] that the TRK/f-sum rule holds in the spin-position basis for self-consistent MBPT. However, there does not exist a proof in the literature that it is satisfied in a lattice system, in which locality gets redefined and the notion of electromagnetic gauge invariance becomes less transparent. In the following, we prove that a response function obtained using equation (12) fulfils the TRK/f-sum rule also in a lattice system. The electromagnetic gauge freedom is described in real space by the transformation (v(r,t), A(r,t)) →(v(r,t)−∂t3(r,t), A(r,t)+∇3(r,t)), where (v(r,t), A(r,t)) and 3(r,t)are the gauge potential and function, respectively. We can extend this transformation to a lattice system by introducing a new gauge function 3i(z), where iis a collective site and spin index. This function is defined on the Keldysh contour of figure 1on which it satisfies the boundary conditions 3i(t0)=3i(t0−iβ). Now, we define the gauge transformation ˆ h(z)→ˆ h3(z)where ˆ h3(z)≡X i j ei(3i(z)−3j(z))hi j (z)ˆc† iˆcj−X i ∂z3i(z)ˆni(A.1) is the gauge transformed one-body part of the Hamiltonian. Moreover, using the equation of motion (4), we can show that G3 i j (z,z0)=ei3i(z)Gi j (z,z0)e−i3j(z0).(A.2) Here, it is found essential that the interaction matrix elements are two-index quantities, since then the exponential phase factors exp(i3i(z)) cancel at each interaction vertex. This cancellation guarantees that (6[G3])i j (z,z0)=ei3i(z)(6[G])i j (z,z0)e−i3j(z0),(A.3) which, in practice, is true only if the one-body Green’s function is calculated self-consistently. The Hamiltonian (A.1) is given, to first order with respect to 3i(z), by h3(z)=h(z)+v(z)+· · · ,(A.4) New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 19 where we defined the one-body potential vi j (z)≡ihi j (z)(3i(z)−3j(z))−δi j ∂z3i(z). (A.5) Consequently, the first-order variation of the one-body Green’s function can be written as δG3 i j (z,z0)=X kl Z γ d¯z Li j,lk (zz0,¯z)vkl (¯z) =iX kl Z γ d¯z(Li j,kl(zz0,¯z)hlk(¯z)−Li j,lk (zz0,¯z)hkl (¯z) −iδkl ∂¯zLi j,lk(zz0,¯z))3l(¯z), (A.6) where Li j,kl(zz0,¯z)is the generalized response function defined in equation (13). The second line above was obtained by using partial integration, as well as the boundary conditions Li j,kl(zz0,t0−iβ) =Li j,kl(zz0,t0)and 3i(t0−iβ) =3i(t0). We can also expand equation (A.2) to first order with respect to 3i(z)and equate the result against the variation (A.6). This leads to a result given by (δikδ(z0,z00)−δjkδ(z,z00))Gi j (z,z0) =X l (hkl (z00)Li j,lk(zz0,z00)−Li j,kl(zz0,z00)hlk (z00)) −i∂z00 Li j,kk (zz0,z00), (A.7) which is nothing but the generalized Ward identity relating the generalized response function to the one-body Green’s function. At the limit z0=z+using equation (14), the Ward identity can be split into three real-time equations, one for each Keldysh component. These equations are given by (δik −δjk)γi j (t)=i(χ> i j,kk(t,t)−χ< i j,kk(t,t)),(A.8a) i∂t0χ≷ i j,kk(t,t0)=X lhkl (t0)χ≷ i j,lk(t,t0)−χ≷ i j,kl(t,t0)hlk(t0),(A.8b) where γi j (t)and χ≷ i j,kl(t,t0)denote the one-body reduced density matrix and the greater/lesser component of the real-time response function, which are defined in equations (6) and (11), respectively. Consequently, the density response function χi j (t,t0)=θ(t−t0)(χ> ii,j j (t,t0)−χ< ii,j j (t,t0))(A.9) satisfies an equation given by i∂t0χi j (t,t0)=δ(t−t0)(χ< ii,j j (t,t0)−χ> ii,j j (t,t0)) +θ(t−t0)X k (hjk(t0)(χ> ii,kj (t0,t) −χ< ii,kj (t0,t)) −(χ> ii,jk(t0,t)−χ> ii,jk(t0,t))hk j (t0)), (A.10) where we used the Ward identity (A.8b). Moreover, evaluating this equation at equal times t0= t−ηwith the limit η→0+taken afterwards and using the symmetry χ> i j,kl(t,t0)=χ< kl,i j (t0,t) (see equation (11)) leads to i∂t0χi j (t,t0)|t0=t−=X k (hjk(t)(χ> kj,ii (t,t)−χ< jk,ii (t,t)) −(χ> kj,ii (t,t)−χ> jk,ii (t,t))hk j (t)). (A.11) New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 20 If the initially unperturbed system is time-translationally invariant, then the density response function depends only on a relative time, or in other words χi j (t,t0)=χi j (t−t0). Then the Ward identity (A.8a) allows us to write the TRK sum, or the f-sum, rule in the time domain as −∂tχi j (t)|t=0+=hji γi j +γji hi j −δi j X k (hjkγkj +γjk hkj ).(A.12) Furthermore, by assuming that the density response function is analytic in the upper half of the complex plane, the TRK/f-sum rule can be written in the frequency domain [20] as −1 πZ∞ −∞ dω ω =m(χi j (ω))=hji γi j +γji hi j −δi j X k (hjkγkj +γjk hkj ).(A.13) Note that in the text, we refer to the right-hand side of these equations as the ground state commutator c1 i j since in an exact theory it can be written [34] as c1 i j = h8N 0|[[ˆni,ˆ H],ˆnj]|8N 0i,(A.14) where |8N 0iis the N-particle ground state, ˆnithe site density operator and ˆ H(t)the manyparticle Hamiltonian. References [1] My¨ oh¨ anen P, Stan A, Stefanucci G and van Leeuwen R 2008 A many-body approach to quantum transport dynamics: initial correlations and memory effects Europhys. Lett. 84 67001 [2] My¨ oh¨ anen P, Stan A, Stefanucci G and van Leeuwen R 2009 Kadanoff–Baym approach to quantum transport through interacting nanoscale systems: from transient to the steady-state regime Phys. Rev. B80 115107 [3] Thygesen K S and Rubio A 2008 Conserving GW scheme for nonequilibrium quantum transport in molecular contacts Phys. Rev. B77 115333 [4] Strange M, Rostgaard C, H¨ akkinen H and Thygesen K S 2011 Self-consistent GW calculations of electronic transport in thiol- and amine-linked molecular junctions Phys. Rev. B83 115108 [5] Baym G and Kadanoff L P 1961 Conservation laws and correlation functions Phys. Rev. 124 287–99 [6] Baym G 1962 Self-consistent approximations in many-body systems Phys. Rev. 127 1391–401 [7] Puig von Friesen M, Verdozzi C and Almbladh C-O 2009 Successes and failures of Kadanoff–Baym dynamics in Hubbard nanoclusters Phys. Rev. Lett. 103 176404 [8] Puig von Friesen M, Verdozzi C and Almbladh C-O 2010 Kadanoff–Baym dynamics of Hubbard nanoclusters: performance of many-body schemes, correlation-induced damping and multiple steady and quasi-steady states Phys. Rev. B82 155108 [9] Strinati G 1988 Application of the Green’s function method to study the optical properties of semiconductors Rev. Nuovo Cimento 11 1–86 [10] Kwong N-H and Bonitz M 2000 Real-time Kadanoff–Baym approach to plasma oscillations in a correlated electron gas Phys. Rev. Lett. 84 1768 [11] Dahlen N E and van Leeuwen R 2007 Solving the Kadanoff–Baym equations for inhomogeneous systems: application to atoms and molecules Phys. Rev. Lett. 98 153004 [12] Starcke J H, Wormit M, Schirmer J and Dreuw A 2006 How much double excitation character do the lowest excited states of linear polyenes have? Chem. Phys. 329 39–49 [13] Romaniello P, Sangalli D, Berger J A, Sottile F, Molinari L G, Reining L and Onida G 2009 Double excitations in finite systems J. Chem. Phys. 130 044108 [14] Sangalli D, Romaniello P, Onida G and Marini A 2011 Double excitations in correlated systems: a many-body approach J. Chem. Phys. 134 034115 New Journal of Physics 14 (2012) 013032 (http://www.njp.org/) 21 [15] Pal G, Pavlyukh Y, Schneider H C and H¨ ubner W 2009 Conserving quasi-particle calculations for small metal clusters Eur. Phys. J. B70 483–96 [16] Dahlen N E and van Leeuwen R 2005 Self-consistent solution of the Dyson equation for atoms and molecules within a conserving approximations J. Chem. Phys. 122 164102 [17] Keldysh L V 1965 Diagram technique for nonequilibrium processes JETP 20 1018 [18] Danielewicz P 1984 Quantum theory of nonequilibrium processes Ann. Phys. 152 239–04 [19] van Leeuwen R, Dahlen N E, Stefanucci G, Almbladh C-O and von Barth U 2006 Introduction to the Keldysh formalism Lect. Notes Phys. 706 33–59 [20] van Leeuwen R and Dahlen N E 2004 Conserving approximations in nonequilibrium Green function and density functional theory The Electron Liquid Paradigm in Condensed Matter Physics, Proc. Int. School of Physics Enrico Fermi vol CLVII (Amsterdam: IOS Press) pp 169–88 [21] Dahlen N E, van Leeuwen R and Stan A 2006 Propagating the Kadanoff–Baym equations for atoms and molecules. J. Phys. Conf. Ser. 35 340–8 [22] Langreth D C 1976 Linear and nonlinear response theory with applications Linear and Nonlinear Electron Transport in Solids ed J T Devreese and V E van Doren (New York and London: Plenum) pp 3–31 [23] Kadanoff L P and Baym G 1962 Quantum Statistical Mechanics (New York: Benjamin) [24] McLachlan A D and Ball M A 1964 Time-dependent Hartree–Fock theory for molecules Rev. Mod. Phys. 36 844–55 [25] Casida M E 1996 Time-dependent density functional response theory of molecular systems: theory, computational methods and functionals Recent Developments and Applications in Density Functional Theory ed J M Seminario (Amsterdam: Elsevier) pp 391–439 [26] Linderberg J and ¨ Ohrn Y 1968 Derivation and analysis of the Pariser–Parr–Pople model J. Chem. Phys. 49 716 [27] Essler F H L, Frahm H, G¨ ohmann F, Kl¨ umper A and Korepin V E 2005 The One-Dimensional Hubbard Model (Cambridge: Cambridge University Press) [28] Saad Y 2003 Iterative Methods for Sparse Linear Systems 2nd edn (Philadelphia, PA: SIAM) [29] Lin H Q and Gubernatis J E 1993 Exact diagonalization methods for quantum systems Comput. Phys. 7400–7 [30] Park T J and Light J C 1986 Unitary quantum time evolution by iterative Lanczos reduction J. Chem. Phys. 85 5870–6 [31] Stan A, Dahlen N E and van Leeuwen R 2009 Levels of self-consistency in the GW approximation J. Chem. Phys. 130 114105 [32] Balzer K, Bauch S and Bonitz M 2010 Time-dependent second Born calculations for model atoms and molecules in strong laser fields Phys. Rev. A82 033427 [33] Stan A, Dahlen N E and van Leeuwen R 2009 Time-propagation of the Kadanoff–Baym equations in inhomogeneous systems J. Chem. Phys. 130 224101 [34] Goodman B and Sj¨ olander A 1973 Application of the third moment to the electric and magnetic response function Phys. Rev. B8200–14 New Journal of Physics 14 (2012) 013032 (http://www.njp.org/)