Kadanoff-Baym approach to double-excitations in finite systems

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 second Born approximations. Our method is equivalent to solving the Bethe-Salpeter equation with a high-level kernel while respecting self-consistently, 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, which is the simplest time-nonlocal approximation, reproduces well most of the additional excitations, which are characterized as double-excitations.


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 . 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 v and gives direct access to the neutral excitation spectra. This quantity satisfies an integral equation known as the 3 Bethe-Salpeter equation (BSE) [9] whose integral kernel is given by the functional derivative δ /δ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 doubleor 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 2 in 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.

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 where . . . denotes an ensemble average with respect to the grand canonical density operator andĉ ( †) 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 z and z along the extended Keldysh contour γ [17] as shown in figure 1; for further details, see [18][19][20][21]. The extended Keldysh contour starting from t 0 and ending at t 0 − 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 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 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 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 is 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 where [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.

5
In particular, all time-local observables are related to the one-body reduced density matrix (1RDM), which is given by where t is 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 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 where ∂ t is proportional to the identity matrix, and tr A ≡ i A ii 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 with respect to an external perturbation v i j (t), and has the general structure of a retarded function given by where the greater and lesser components are defined, respectively, in terms of the 1RDM operatorγ H,i j (t) ≡ĉ † H, j (t)ĉ H,i (t), as 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 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]. The Bethe-Salpeter equation for the generalized response function whose particular time ordering L i j,kl (zz + , z ), 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.

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 and in terms of which the greater and lesser components of the real-time response function, defined in equation (11), are given by 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 where we defined the four-point Bethe-Salpeter kernel as the functional derivative Figure 2 contains 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 , corresponds to a solution of the Bethe-Salpeter equation with the kernel δ /δ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). Figure 3. The Hartree-Fock and second Born self-energies ( HF/2B ) and Bethe-Salpeter kernels (K HF/2B = δ HF/2B /δG HF/2B ). Each self-energy diagram shown relates to at least one kernel diagram, for example, bubble (B) gives rise to diagrams B 1 − B 3 . Wiggly lines denote the Coulomb interaction, and double lines denote dressed Green's functions.

Approximate Bethe-Salpeter kernels
The time propagation of a Green's function with an approximate self-energy automatically leads to a retarded response function related to an approximate kernel δ /δ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 n Green's function lines gives rise to n inequivalent 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 L i j,kl (zz + , z ). 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 B 1−3 and X 1−3 can be deciphered similarly. Diagram B 1 describes a partially screened interaction between a particle and a hole, while diagrams B 2 -B 3 reflect the  . Some diagrams containing multiple simultaneous particle-hole pairs, whose origin lies in either self-energy insertions (a) or kernel diagrams (b). A simple line represents a Hartree-Fock Green's function in terms of which the top and bottom rows describe two-and three-particle excitations, respectively. The direction of time is from the left to the right.
exchange symmetry or describe processes in which the original particle-hole pair is annihilated and another emerges in the event. Diagrams X 1−3 relate only to partial exchange processes in which either an entering particle (X 1,2 ) or hole (X 1,3 ) gets annihilated and another particle or hole created during the process exits.
We will now discuss the effect of the temporal structure of an approximation to the excitation spectrum which is given by the pole structure of the response function. We know that a time-local kernel only reproduces poles whose locations are given by differences in the renormalized, single-particle excited and ground state energies [13,14,24,25]. The excitations corresponding to these poles are then called one-particle excitations, while all remaining ones we refer to as many-particle excitations. However, we can also view these one-to many-particle excitations in terms of some single-particle basis which has an appropriate energy label. In this case, excitations can have one-to many-particle character depending on their representation in this basis. In a wave-function-based theory, excitations are characterized with the help of the basis expansion of the many-particle excited and ground states. On the other hand, in many-body perturbation theory, we examine the diagrammatic structure of a response function in terms of Green's functions represented in this single-particle basis.
In figure 4 we show some low-order diagrams which are obtained by expanding 2B Green's functions in terms of HF Green's functions. These diagrams represent time-non-local processes that contain multiple, simultaneous particle-hole pairs, or in other words describe many-particle excitations in the HF basis. Such additional excitations are clearly not present within the HF approximation and can, in general, arise either from self-energy insertions or directly from the kernel as shown in figures 5(a) and (b), respectively. Moreover, self-consistency ensures that some classes of these diagrams representing many-particle excitations are taken into account to all orders of perturbation theory. If some of these diagrams are dominant in the perturbation expansion, we can expect many-particle excitations to appear in the excitation spectra.

Model and test systems
We consider N -electron lattice systems which are based on the Pariser-Parr-Pople (PPP) model [26]. This model is commonly used to study properties of organic molecules, and in 9 (a) (b) Figure 5. The red sites in the spatial profiles of our test systems, shown in (a), denote the target site for a perturbation, while d ≡ d i j = 1 is the distance between nearest neighbors. In (b), we show for the interaction strengths U = 1.0, 2.0 the HF single-particle energies for the spin-restricted ground state in which the lowest energy levels are occupied. Note that the energy spectra of the ring coincide with the non-interacting spectra up to a constant because of the discrete rotational symmetry.
the limit of short interactions, it reduces to the Hubbard model. The Hubbard model [27] is a standard model to describe strong correlations in condensed matter physics. Our Hamiltonian is given byĤ whereĉ ( †) iα annihilates (creates) an electron of spin α at site i, whilen i ≡ĉ † iσĉ iσ andn i ≡ αn iα are site density operators. The one-body Hamiltonian is given by where the hopping elements δ i, j are equal to 1 for the nearest-neighbor sites on the lattice and are zero otherwise. The symbol v i (t) = vδ 1,i δ(t) denotes a time-dependent potential. The two-body, electron-electron interaction is described with local and non-local interaction matrix elements denoted by U , known as Hubbard-U , and w i j , such that w ii ≡ 0, respectively. These interaction matrix elements and the positive background potential Z are chosen to address two kinds of systems itemized below: (i) The Hubbard model is specified by choosing the off-diagonal interaction matrix elements, as well as the positive background potential, as zero (w i j = 0, Z = 0). The model represents a system with a short range or more precisely a contact two-body interaction. (ii) The PPP model, on the other hand, is introduced to address the effects of long-range interactions, which are described by the Coulomb-type matrix elements w i j = U/2d i j , where d is a parameter describing the distance between sites i and j. We set the distance between nearest neighbors on a lattice, as well as the positive background potential, to one (d i, j = 1, Z = 1).
We have, in particular, prepared lattice systems of two main types: linear chains and uniform rings, and performed calculations by varying the interaction strengths and the number of electrons, but restricted ourselves to spin-compensated systems. In the next section, we summarize some of our results by considering a six-site, six-electron (half-filled) chain with long-range interactions, as well as a six-site, two-electron (one-sixth-filled) ring with shortrange interactions. Some essential features of both test systems are represented in figure 5. Note that we do not aim to model a realistic physical system, but instead we have chosen these energetically simple systems for purely benchmark purposes.

Results: addition, removal and excitation spectra
We have chosen two finite lattice systems, described in section 3, to test the performance of the many-body approximations obtained from the HF and 2B self-energies. The observables used as a measure of the quality of these approximations are the ground state energy, and the non-neutral and neutral excitation spectra calculated using the spectral and retarded response functions, respectively. The quality of our many-body schemes is assessed by comparing all approximate quantities against exact results.

Numerical aspects
All exact quantities are obtained by exact diagonalization (ED), or in other words calculated using the many-particle eigenstates of the system which are obtained by solving the time-(in)dependent Schrödinger equation. In the spin-compensated subspace, the size of the manyparticle basis is given by N S N /2 2 , where N S and N are the number of sites and electrons, respectively. Our numerical implementation of ED in these finite lattice systems relies on the well-known Lanczos routine [28][29][30]. Approximate quantities, on the other hand, are obtained by calculating the equilibrium Green's function [16,31] which yields the ground state energy and from which the non-equilibrium Green's function is obtained by time propagation. These calculations are performed with a parallelized version [32] of the algorithm described in [21,33].
Some additional remarks on the calculation of spectral properties are in order. These properties are, in practice, computed either by direct evaluation of addition, removal and excitation energies (only in ED) or by finite-length time propagation of a many-particle ground state (ED) and an equilibrium Green's function (MBPT). The latter approach leads to timedomain spectral and response functions which are transformed into their frequency-domain counterparts. As a side-effect of the finite propagation length, oscillatory noise-like features, as well as broadening of peaks, appear in these frequency-domain spectra. Such defects can be reduced by increasing the propagation length denoted by T and reached with a time step . In figures, results related to the direct evaluation scheme will be denoted by ED-D, while time-propagated results are denoted by ED, HF or 2B.

Addition and removal energies
The spectral function relates to the so-called quasi-particle properties, in particular the frequency-domain structure contains information about particle addition and removal energies, and the diagonal elements represent the local density of states. The Lehmann representation of the ground state spectral function is given by where g N +1 n,iα ≡ N 0 |ĉ iα | N +1 n − E N 0 are the addition and removal energies. The symbols | N k and E N k denote the kth eigenstate and eigenenergy of an N -particle system, respectively.
The spectral functions of the half-filled chain with different interaction strengths (U = 1.0, 2.0) are shown in figure 6. Both many-body approximations capture well the lowest addition and removal energies that clearly relate to the lowest unoccupied and highest occupied single-particle states. This happens irrespective of the interaction strength. Moreover, the approximations open up the band gap adequately and reproduce correctly the particle-hole symmetry, which although being in general broken by the long-range interaction has been restored by the choice Z = 1 for the positive background. Differences between approximate solutions, as well between approximate and exact spectra, appear at higher energies; see the addition energies/removal energies A2-A5/R2-R5. These energies can be only partly related to renormalized single-particle states, as correlation effects give rise to entirely new structure. Our results indicate that the 2B approximation reproduces some of this additional structure present in the exact spectral function, while the HF method, merely capturing the single-particle energy spectra (Koopmans' theorem), fails to do so.
The addition and removal energy spectra of the one-sixth-filled, six-site (two electrons) ring are given in figure 7. When U = 1.0 both approximations capture well the lowest addition energies and removal energy, which relate again to the lowest unoccupied and highest occupied states. Additionally, they reproduce correctly the broken particle-hole symmetry, irrespectively of the interaction strength. As the interaction strength is increased, the HF approximation underestimates both the removal energy and overestimates all addition energies   which, altogether, signals a failure to reproduce both the two-electron ground state and the lowest three-electron states correctly. The ground state energies given in table 1 verify these speculations, as well as the fact that both of our approximations describe the one-particle system correctly and therefore the removal energies describe solely the quality of the ground state. Although both approximations reproduce the addition energies (A1, A3) with the highest intensities, our correlated approximation gives more accurate results. Moreover, only the 2B approximation is able to quantitatively account for the additional, low-intensity structure which consists of a pole between high-intensity ones (A2) and two higher-lying poles (A4, A5). These energies correspond to the addition of an electron to the highest unoccupied single-particle state (A4), as well as to more complicated correlation-induced states (A2, A5). Overall, we have observed that the HF approximation gives quantitatively incorrect highenergy spectra, which is simply due to the fact that there are fewer single-particle energies than eigenenergies of the true interacting system. The 2B approximation, on the other hand, capable of accounting for the so-called correlation-induced quasi-particles, does reproduce some of the warranted additional structure, and additionally captures correctly shifting intensities as interaction strength is increased; see poles A2/R2 of figure 6 or A2 of figure 7. In the next section, we will address the question of how well these fundamental differences affect the excitation structure evaluated at a higher δ /δG level.

Excitation energy spectra
The neutral excitation energies can be accessed through the density(-density) response function, which is given by where the retarded response function is defined in equation (9). In the following, we consider time-translationally invariant systems in which the response function depends only on a relative time; that is, χ iα, jβ (t, t ) = χ iα, jβ (t − t ). The connection to the excitation energies becomes explicit in the Lehmann representation, in which the imaginary part of the ground state density response function can be written as where f n,iα ≡ N 0 |n iα | N n is the oscillator strength and n ≡ E N n − E N 0 the neutral excitation energies. Such an equation is valid as long as the eigenstates of the N -particle system, and hence the oscillator strengths, can be chosen real valued, which is the case for our time-reversal invariant systems.
However, in practice, we calculate the ground state density response function χ i j (t) ≡ αβ χ iα, jβ (t) by means of time propagation. First, we perturb the system with a spatially and temporally local potential v i (t) = v j δ i j δ(t), where j denotes the target site for the perturbation, v is the strength of the perturbation and δ(t) is a Dirac delta function. Then, we use equation (12), which states that the density response function is given by where n 0 i denotes the ground state site density, and n i (t) = α γ iα,iα (t) is the site density of the perturbed system evaluated at time t. Although nonlinear effects are always present, they are easily reduced by ensuring that the magnitude of the perturbation lies well in the linear response region. In practice, this is verified by doubling the perturbation and checking that also the response doubles to a sufficient accuracy.
Before presenting neutral excitation spectra, we will introduce some auxiliary quantities that are useful for the discussion of the results.
(i) The ground state energy can be used to judge whether a possible deviation from an exact excitation energy results from a failure to describe the ground or excited state correctly. Although we cannot access the ground state energy at the δ /δG level, a first-order approximation denoted by E GM is provided by the GM functional (8) evaluated with the 2B self-energy. (ii) We prove in the appendix that any diagrammatically well-defined, self-consistent manybody approximation satisfies the f-sum rule; that is, the density response function satisfies equations (A.12) and (A.13). In our case, this simply means that where h i j and γ i j ≡ αβ γ 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 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 bŷ whered ( †) 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 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 where Tr is the trace over a complete set of many-particle states, and a properly normalized density operator which projects any state to a degenerate subspace, with eigenvalue E, of the many-body eigenstates. The oscillator strength f k,iα is defined below equation (21). In practice, we calculate the differences where E 0 is 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 8 for 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  Figure 9. The relative weighted occupation numbers, N k (E) ≡ N HF 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 E GM -E 0 ≈ 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 9 for 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 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.

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 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.
where we defined the one-body potential v i j (z) ≡ ih i j (z)( i (z) − j (z)) − δ i j ∂ z i (z). (A.5) Consequently, the first-order variation of the one-body Green's function can be written as where L i j,kl (zz ,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 L i j,kl (zz , t 0 − iβ) = L i j,kl (zz , t 0 ) and i (t 0 − iβ) = i (t 0 ). We can also expand equation (A.2) to first order with respect to i (z) and equate the result against the variation (A.6). This leads to a result given by (δ ik δ(z , z ) − δ jk δ(z, z ))G i j (z, z ) = l (h kl (z )L i j,lk (zz , z ) − L i j,kl (zz , z )h lk (z )) − i∂ z L i j,kk (zz , z ), (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 z = 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 where γ i j (t) and χ ≷ i j,kl (t, t ) 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 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, t ) = χ i j (t − t ). 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 + = h ji γ i j + γ ji h i j − δ i j k (h jk γ k j + γ jk h k j ). (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 Note that in the text, we refer to the right-hand side of these equations as the ground state commutator c 1 i j since in an exact theory it can be written [34] as 14) where | N 0 is the N -particle ground state,n i the site density operator andĤ (t) the manyparticle Hamiltonian.