The G1--G2 Scheme: Dramatic Acceleration of Nonequilibrium Green Functions Simulations Within the Hartree--Fock-GKBA

The time evolution in quantum many-body systems after external excitations is attracting high interest in many fields. The theoretical modeling of these processes is challenging, and the only rigorous quantum-dynamics approach that can treat correlated fermions in two and three dimensions is nonequilibrium Green functions (NEGF). However, NEGF simulations are computationally expensive due to their $T^3$-scaling with the simulation duration $T$. Recently, $T^2$-scaling was achieved with the generalized Kadanoff--Baym ansatz (GKBA), for the second-order Born (SOA) selfenergy, which has substantially extended the scope of NEGF simulations. In a recent Letter [Schl\"unzen \textit{et al.}, Phys. Rev. Lett. \textbf{124}, 076601 (2020)] we demonstrated that GKBA-NEGF simulations can be efficiently mapped onto coupled time-local equations for the single-particle and two-particle Green functions on the time diagonal, hence the method has been called G1--G2 scheme. This allows one to perform the same simulations with order $T^1$-scaling, both for SOA and $GW$ selfenergies giving rise to a dramatic speedup. Here we present more details on the G1--G2 scheme, including derivations of the basic equations including results for a general basis, for Hubbard systems and for jellium. Also, we demonstrate how to incorporate initial correlations into the G1--G2 scheme. Further, the derivations are extended to a broader class of selfenergies, including the $T$ matrix in the particle--particle and particle--hole channels, and the dynamically screened-ladder approximation. Finally, we demonstrate that, for all selfenergies, the CPU time scaling of the G1--G2 scheme with the basis dimension, $N_b$, can be improved compared to our first report: the overhead compared to the original GKBA, is not more than an additional factor $N_b$.


I. INTRODUCTION
Nonequibrium Green functions (NEGF) [1][2][3], have proven highly successful in simulations of the dynamics of correlated many-body systems. This is due to a number of attractive properties that include conservation laws and the existence of systematic approximations schemes that are based on Feynman diagrams. Moreover, NEGF allow for a rigorous derivation of quantum kinetic equations and for their systematic improvement; for recent overviews, see the text books [4][5][6].
While early computational applications focused on spatially homogeneous systems such as nuclear matter [7,8], optically excited semiconductors [4,9], and dense plasmas [10,11], during the recent 15 years the scope of applications has substantially broadened. This includes the excitation and ionization dynamics of small atoms and molecules [12][13][14], the correlated-electron dynamics in the Hubbard model [15][16][17], the dynamics of fermionic atoms [18,19], and the stopping of ions in correlated materials [20][21][22]. This success was caused, among others, by progress in the numerical solution of the basic equations of NEGF-the Keldysh-Kadanoff-Baym equations [12,[23][24][25]. Furthermore, improved time propagation and integration schemes have allowed to increase the efficiency and accuracy of the simulations [26,27]. Moreover, the implementation of more advanced selfenergies, such as the T -matrix selfenergy, have allowed to increase the accuracy and predictive capability; for a recent review, see Ref. [28]. In particular, very good agreement with cold-atom experiments [18] and with ab initio densitymatrix-renormalization-group (DMRG) simulations were reported [27]. A particular advantage of NEGF simulation is that they are well capable to treat electronic correlations, in contrast to density-functional theory (DFT), and that they are neither restricted to 1D systems, such as DMRG, nor to short times, such as continuous-time quantum Monte Carlo [29].
The main disadvantage of NEGF is their high numerical effort. The majority of many-body methods, including time-dependent DFT (TDDFT), Boltzmann-type quantum kinetic equations, hydrodynamics or semiclassical molecular dynamic-and even the exact solution of the time-dependent Schrdinger equation-require a simulation time that grows linearly with the physical time. In contrast, for NEGF, the propagation in the two-time plane, together with the memory integration in the scattering contributions, gives rise to a N 3 t -scaling, where N t is the propagation time (number of time steps). A substantial acceleration is possible when the generalized Kadanoff-Baym ansatz (GKBA) is applied [30] which restricts the propagation to a time-stepping along the time diagonal. If combined with Hartree-Fock propagators (HF-GKBA) [31][32][33] the CPU time scaling can be reduced to N 2 t , which has given rise to a drastic increase of the number of HF-GKBA simulations in recent years, e.g. Refs. [17,27,[34][35][36][37][38]. However, this improved scaling is achieved only for the simplest selfenergy-for the second-Born approximation (SOA). If the HF-GKBA is applied to improved selfenergies, such as the T -matrix selfenergy [19,27], which is required for strongly corre-lated systems [18], or the GW selfenergy [28] which is required to capture dynamical screening effects, the CPU time scaling is again increased to N 3 t .
In a recent Letter we reported a breakthrough for NEGF simulations within the HF-GKBA scheme: we demonstrated that time-linear scaling, i.e. a CPU time that is of order N 1 t , can be achieved if the equations of motion are reformulated, without any approximations. The alternative approach solves the time-local equations for the time-diagonal single-and two-particle Green functions and was called G1-G2 scheme [39]. While the equivalence of the HF-GKBA to time-local equations was pointed out before [5,40], a comparison of the numerical behavior of both approaches was performed only in Ref. [39]. There we predicted N 1 t -scaling for SOA and GW selfenergies and any type of single-particle basis. The scaling was demonstrated for small Hubbard clusters which turned out to be the most unfavorable case because the CPU time of the G1-G2 scheme was found to grow by a factor N 2 b faster than for the standard HF-GKBA approach.
In this article we present extensive additional results for the G1-G2 scheme. First, we present all necessary details for the derivation of the equation of motion for the time-diagonal two-particle Green function. The results are derived for a general basis, for the Hubbard model and for jellium. Second, we discuss how initial correlations can be incorporated. Third, we extend the analysis to other selfenergies: the T -matrix approximation in the particle-particle (TPP) and particle-hole (TPH) channels and the dynamically-screened-ladder (DSL) approximation. Fourth, numerical results are demonstrated for all selfenergy approximations which clearly confirm the N 1 tscaling. Finally, we re-evaluate the N b -dependence of the CPU time and report an additional optimization that reduces the overhead of the new scheme from N 2 b to only N 1 b , for the Hubbard model, for all selfenergies.
This paper is organized as follows: In Sec. II we summarize the main required formulas of NEGF theory and the properties of the two-particle Green function. In Sec. III we present the basic formulas for the G1-G2 scheme for the case of SOA selfenergy-separately, for a general basis, the Hubbard basis and for jellium. The same analysis is then extended to GW and T -matrix selfenergies in Secs. IV and V, and to the screened-ladder approximation in Sec. VI. Finally, the analysis of the scaling behavior with N t and N b for all selfenergies and numerical results are presented in Sec. VII.

II. THEORETICAL FRAMEWORK A. Keldysh-Kadanoff-Baym Equations and two-particle Green function
We consider a nonequilibrium quantum many-particle system with the generic Hamiltonian containing a single-particle contribution h (0) and a pair interaction w. The matrix elements are computed with an orthonormal system of single-particle orbitals |i . The creation (ĉ † i ) and annihilation (ĉ i ) operators of particles in state |i define the one-body nonequilibrium Green function (correlation function) for contour-time arguments z on the Keldysh contour C [28], Here, T C is the time-ordering operator on the contour, and the averaging is performed with the correlated unperturbed density operator of the system. The equations of motion (EOMs) for the NEGF are the Keldysh-Kadanoff- (the times z ± := z ± , 1 are slightly shifted to disambiguate time ordering) that couple to the two-particle Green function which contains a mean-field (Hartree-Fock) and a correlation contribution Our scheme involves the special case of two-particle functions that depend either on one or two times and their real-time components [28], that we define as follows The time-diagonal two-particle Green function, G(t), defined by Eq. (4), is the central quantity of the G1-G2 scheme. In general, and for all considered selfenergy approximations in this work, it obeys the following (pair-) exchange symmetries, B. Time-diagonal KBE Following Eqs. (1) and (2) the EOM for the less component of the NEGF on the real-time diagonal, G ≷ ij (t) := G ≷ ij (t, t), can be rewritten as [42] i d dt where the right-hand side contains the collision integral Here, t 0 marks the time, at which the system's initial state is prepared; the treatment of initial correlations is discussed in Sec. III F. On the time diagonal the less component of the NEGF can be written as where n ij is the single-particle density matrix. In Eq. (7) the mean-field part of the two-particle Green function, cf. Eq. (3), is included in an effective single-particle Hartree-Fock Hamiltonian which is defined as [28] h HF where we introduced the (anti-)symmetrized fourdimensional interaction matrix elements The interaction tensor obeys the same symmetries as G(t) [cf. Eqs. (5) and (6)], which also leads to The selfenergy Σ in the collision integral (8) contains only the remaining correlation contribution of the two-particle Green function.

III. SECOND-ORDER BORN SELFENERGY
In the following we introduce the G1-G2 scheme for the simplest case of choosing the selfenergy in the second-Born approximation [42], With that, the collision integral of the time-diagonal equation (8) transforms into: where we presented several equivalent formulations that will be used below. At this point, it is possible to identify where, in the transformation (*), the symmetry of Eq. (13) is used and two summation indices are switched, r ↔ s.
A. G within the GKBA The G1-G2 scheme is a reformulation of the ordinary solution of the time-diagonal KBE in the HF-GKBA. When applying the GKBA the time off-diagonal elements of the less and greater NEGF are reconstructed from the time diagonal value via [42] We now rewrite this reconstruction, taking into account that in the collision integral (8) only G > (t ≥t) and G < (t ≤ t) appear. Replacing them with the GKBA and taking into account the special case of equal times results in where we introduce the modified propagator (timeevolution operator) which on the time diagonal reduces to [42], By inserting the GKBA, Eqs. (16) and (17), for the Green functions into Eq. (14), we find (cf. Appendix A 1) where we now introduce short notations for the twoparticle propagator U (2) and define the occupation factors Φ ≷ , An even more compact notation can be achieved by introducing the two-particle source term which results in The function Ψ ± has the meaning of pair correlations produced in the system via two-particle scattering per unit time.

B. Time-linear integral solution for G
When applying the HF-GKBA the retarded and advanced propagators, G R (t, t ) and G A (t, t ), are described in HF approximation and, thus, obey the group property [36] for t >t > t : As we show in Appendix A 2 it directly follows that the two-particle propagator also possesses the group property, ijpq (t,t)U rskl (t, T + ∆) . Applying the group property of the two-particle propagator, Eq. (24), leads to xyjm (T, T + ∆) , where we identify the two-particle Green function at time T , The above expression only contains a time integral of fixed length ∆. Thus, provided that the solution G (T ) is known, the propagation to T +∆ can be done in a constant amount of time, independent of T . While Eq. (25), in principle, provides the basis for a time-linear propagation scheme, its integral form proves to be unfavorable for numerical implementation. Therefore, in the following, an alternative approach (G1-G2 scheme [39]), that is based on the solution of a (time-local) differential equation, is derived which will be analyzed throughout this paper.
C. Time-linear differential solution for G : SOA-G1-G2 equations for a general basis Here we consider a general single-particle basis where spin degrees of freedom are included in the basis index. Below we will separately consider the special cases of a Hubbard basis and the jellium model for electrons where the two spin projections will be indicated explicitly. In order to find the differential equation for G , the EOMs for the retarded/advanced Green functions in HF-GKBA along both time-directions and the diagonal are repeated [42]: For the two-particle propagators similar Schrdinger-type EOMs hold as shown in Appendix A 3, where we define the two-particle Hartree-Fock Hamiltonian as the sum of two single-particle parts: With that we now compute the time derivative of the time-diagonal two-particle Green function within the HF-GKBA, G , which contains two parts, .
The first contribution ( ) originates from the integration boundaries, where the latter equation holds due to the identity [cf.
Eqs. (19) and (20)] The second contribution to the derivative results from the time dependence of the integrand, i.e. of U (2) , and, using the results from Eqs. (27) and (28), we obtain where we identify G again, to get With that, the full derivative of the time-diagonal twoparticle Green function is obtained by adding up the results of Eqs. (30) and (31), We now summarize the equations of the G1-G2 scheme for second order Born selfenergies, for a general basis. The scheme consists of the equation for the time-diagonal element of the single-particle Green function, cf. Eq. (7), coupled to Eq. (32)-the EOM of the time-diagonal element of the two-particle Green function. Equations (32), (33), and (34) constitute a closed system of time-local differential equations, for which the computational effort for a numerical implementation scales linearly with time. This was achieved by eliminating the non-Markovian (memory) structure of the collision integral. All transformations so far are exact and reproduce the standard HF-GKBA result, as was demonstrated in Ref. [39]. The linear scaling with N t , as opposed to the quadratic scaling of the standard HF-GKBA in SOA, is the basis for a potentially dramatic speedup of NEGF simulations. The price to pay is the need to compute the entire matrix of the time-diagonal two-particle Green function the effort for which only depends on the basis dimension N b . This will be analyzed in detail in Sec. VII. In similar manner as for the SOA selfenergy, a timelocal equation for G corresponding to more advanced selfenergies can be derived for which the speedup of the G1-G2 scheme is even larger. This will be demonstrated in the subsequent sections. But before that, we consider the G1-G2 scheme in SOA for two important special cases of basis sets-the Hubbard basis and the spatially uniform jellium model (plane-wave basis).

D. SOA-G1-G2 equations for the Hubbard model
The Hubbard model [43] is among the fundamental models in condensed matter physics, in particular, for the analysis of strong electronic correlations. More recently it has been widely used to study the behavior of fermionic and bosonic atoms in optical lattices [44] and, in particular time-dependent correlation phenomena, see, e.g, Refs. [16,17,45,46]. For the Fermi-Hubbard model, the general pair-interaction matrix element becomes with the on-site interaction U and the spin projection labeled by greek indices. The kinetic energy matrix is replaced by a hopping Hamiltonian, which includes hopping processes between nearestneighbor sites i, j with amplitude J. Thus, the total Hamiltonian is given bŷ Extensions to more complicated models, going beyond the nearest neighbor single-band case are straightforward but will not be considered here. The time-diagonal EOM for the single-particle Green function, Eq. (7), takes the following form (from here we give all Hubbard equations for the spin-up component; the spin-down equations follow from the replacement ↑ ↔ ↓.) where for electrons there exist two collision integrals, I ↑ and I ↓ , that enter the single-particle EOMs.
The equation for the time-diagonal two-particle Green function, Eq. (32), now reads where and The Eqs. (37) and (39)  As the second example we consider the jellium Hamiltonian [47], with the Coulomb matrix element v |q| = 4πe 2 |q| 2 . This model is of relevance for the electron gas in metals [48,49], for electron-hole plasmas in semiconductors [4], and for dense quantum plasmas, e.g. [5,10], as well as for model development [47,49].
The matrix element of the pair interaction in a plane wave basis is where v q denotes the spatial Fourier transform of the pair potential, and the delta function arises from momentum conservation (spatial homogeneity).
The EOM for the single-particle Green function, Eq. (7), is now where we defined and the equation for the time-diagonal two-particle Green function becomes where h HF,α This result agrees with the one derived in Refs. [5,32].

F. Initial pair correlations in the G1-G2 scheme
We conclude this section by considering the question of initial values in the G1-G2 scheme. Obviously, the solution of the differential equations (33), for G < (t), and Eq. (32), for G(t), are defined only up to arbitrary constants which we can fix by choosing the initial values, G 0,< = G < (t = t 0 ) and G 0 = G(t = t 0 ). Recalling the definitions (9) and (4), the former is related to the initial value of the single-particle density matrix, and the latter to the initial value of the correlated part of the two-particle density matrix, i.e., pair correlations existing in the system at the initial time t = t 0 . While, mathematically, any initial value is compatible with the differential equation, physical considerations do impose restrictions, as was discussed e.g. in Refs. [21,50]. The result can be summarized as follows: only such pair correlations are physically relevant that can be produced by a dynamic evolution of the form starting from an uncorrelated system att → −∞. We underline that the treatment of initial values in the G1-G2 scheme is not restricted to the second-Born approximation but can be generalized to more sophisticated selfenergies.
In the context of NEGF theory and the GKBA, the question of initial correlations has been extensively discussed before, see, e.g., Refs. [6,7,11,51], for more recent investigations, see Refs. [21,36,52]. Going back to our starting point-the integral representation of G, cf. e.g. Eqs. (14) and (21)-it is clear that these expressions vanish, in the limit t → t 0 , i.e., these expressions are valid only for the case of an initially uncorrelated system. These integral solutions are readily extended to the case of arbitrary initial correlations [5]: in that case, the previous solution, Eq. (21), has to be supplemented by a homogeneous solution of the differential equation for G, which we denote G IC , which recovers the structure of Eq. (25). Both terms (47) can be combined into a single expression according to While Eq. (48) in its presented form holds for the second-Born approximation, this functional form is generally valid. The main difference, for more complicated selfenergies, is the explicit form of the two-particle propagators. For the additional approximations considered in this work [GW (Sec. IV), T matrix (Sec. V)] the respective expressions are presented in Appendix C.

IV. GW SELFENERGY
The static second-Born approximation that was considered above neglects screening effects and the dynamics of screening. These effects are captured by the GW approximation for which the selfenergy is given by, Here, W is the dynamically screened interaction, which can be expressed in terms of the bare interaction and the inverse dielectric function, which allows us to transform the selfenergy (49) into, The collision integral of the time-diagonal equation then becomes, Recalling the definition (8), we identify the time-diagonal element of the two-particle Green function in GW approximation, By construction, the screened interaction tensor obeys the following symmetry [cf. Eq. (11)], From Hedin's equations [53] we derive the following relation for the dynamically screened interaction W from which we subtract the singular part, By comparison with Eq. (50) and using the symmetry of Eq. (51) one can identify a recursive equation for ε −1 , The time-diagonal equation for the inverse dielectric function can be further simplified, A. GW approximation within the HF-GKBA We now apply the HF-GKBA [cf. Eqs. (15) - (18)] and obtain the following expressions for G , as well as for ε −1 , where U is given by Eqs. (A2) and (A3). By using the symmetry relation of Eq. (51) we easily find an expression for the time derivative of the off-diagonal inverse dielectric function, where we introduced the modified two-particle Hartree-Fock Hamiltonian which matches the index structure of the effective quasi-Hamiltonian, defined as Combining these Hamiltonians into a single one, we observe that the inverse dielectric function, within the GW -HF-GKBA, obeys a time-dependent two-particle Schrdinger equation, with the Hamiltonian (56), that is equivalent to the rather complicated integral equation (54).
In the following, we demonstrate that, for the GW -HF-GKBA, again, a time-local G1-G2 scheme can be derived which retains time-linear scaling [39].

B. GW -G1-G2 equations for a general basis
To derive the G1-G2 scheme, we compute the time derivative of G , yielding, where the first contribution, which originates from the derivative of the integration boundaries, is given by Here, the two-particle source term is defined as The second contribution to Eq. (58), resulting from the time derivative of ε −1 , is given by whereas the third contribution to Eq. (58), which stems from the derivative of the propagators, is Finally, the three contributions to the derivative of G are combined to reveal where h ε (t) was defined in Eq. (56). With this we have obtained the equations of the G1-G2 scheme for the GW approximation. For h ε,corr (t) ≡ 0, we recover the equations from the SOA, cf. Eq. (32), since the remaining Hamiltonian contribution can be expressed as a commutator. Equation (59) is the most compact formulation that visualizes the intrinsic structure of G in the GW approximation.
For practical use, it is convenient to separate the correlation contributions from the mean-field terms via the introduction of an additional quantity: where polarization effects are included in Equation (59) agrees with the polarization approximation of density-matrix theory, cf. Refs. [5,54]. In the Markov limit this leads to the quantum generalization of the Balescu-Lenard kinetic equation [55][56][57].
Here, we have employed the standard definition of GW in NEGF theory, which is widely used in literature (see, e.g., Refs. [15,28,58]), in which the screened interaction [(52)] does not include exchange terms. The generalization to also describe exchange processes is, however, straightforwardly carried out. For the G1-G2 scheme, this is achieved by simply replacing Ψ ijkl (t) by Ψ ± ijkl (t) in Eqs. (59) and (60).
Again we have succeeded to eliminate all time integrations which means that Eq. (59) can be solved with an effort that is first order in N t . Note that the conventional HF-GKBA scheme with GW selfenergy scales as N 3 t indicating a huge advantage of the G1-G2 formulation [39]. More computational details will be given below, in Sec. VII.   (35). With that, the equations of motion (60) become, where we introduced the polarization terms, Notice that there are two separate spin combinations (four when considering ↑ ↔ ↓) for the two-particle Green function that enter Eqs. (62) and (63). Due to the crosscoupling in the two polarization terms, they cannot be solved independently [28,59]. Numerical results for the GW -G1-G2 scheme are presented in Sec. VII.
With that, the equation (60) for the time-diagonal twoparticle Green function [recall the definition (45)] be-comes, with the momentum representation of the polarization term, given by As we will discuss in Sec. VII, the GW equations for jellium can be solved particularly efficiently.

V. T -MATRIX SELFENERGIES
We next turn to the case of strong coupling where the second-Born approximation is not applicable. It is well known that the entire Born series can be summed up, giving rise to the T -matrix (or binary-collision or ladder) approximation. Here we first consider the case of a static pair interaction. The extension to a dynamically screened T matrix will be considered in Sec. VI. We start by considering, in Sec. V A, the T matrix in the particleparticle channel after which we analyze, in Sec. V B, the T matrix in the particle-hole channel.
A. T matrix in the particle-particle channel For the particle-particle T matrix, the selfenergy has the form [3,19], Here, the T matrix is expressed as which allows us to rewrite the selfenergy (67): In Eqs. (68) and (69) the quantity Ω pp is the nonequilibrium generalization of the Møller operator from scattering theory [60,61]. The collision integral (8) of the time-diagonal equation then becomes, which results in the following expression for the timediagonal element of the two-particle Green function, By construction, the T matrix obeys the following symmetry [cf. Eq. (12)], The T matrix sums up the particle-particle collisions via the recursive equation (nonequilibrium Lippmann-Schwinger equation; compared to the standard definition of the T matrix, here the singular part has been subtracted [19,28]), Following this and using the symmetries of Eqs. (12) and (70) the relation for the Møller operator is readily derived, The time-diagonal equation for Ω pp can be further simplified,

T pp approximation within the HF-GKBA
We now apply the HF-GKBA [cf. Eqs. (15) - (18)] and find the following expressions for G , as well as for Ω pp , As for the case of GW selfenergies, here we introduced two quasi-Hamiltonians, Combining these Hamiltonians again into a single one, the equation (72) for the Møller operator is transformed into a time-dependent two-particle Schrdinger equation, This equation is analogous to the Schrdinger equation for the inverse dielectric function, Eq. (57), the main difference being the modified Hamiltonian (74).

T pp -G1-G2 equations for a general basis
To derive the G1-G2 scheme for the particle-particle T matrix, we have to take the derivative of G, yielding, The derivative of the integration boundaries results in, while the time derivative of the Møller operator yields, The last contribution originates from the derivative of the two-particle propagator, Combining the three contributions to the derivative of G reveals where h Ω pp (t) was introduced in Eq. (74). This is the central equation for the G1-G2 scheme in T -matrix approximation for the particle-particle channel [5,61]. Compared to the equation of motion for G in second Born approximation, Eq. (32), this equation contains, in addition, the particle-particle ladder terms which are generated by the quasi-Hamiltonian of Eq. (73). Again, for practical use, it is convenient to separate the correlation contributions from the mean-field terms via the introduction of an additional quantity: where the particle-particle ladder term is defined by Without the Λ-terms we exactly recover the equation of motion for G in second-order Born approximation. Inclusion of the Λ-terms, on the other hand, allows one to take into account multiple scattering and large-angle scattering effects that are important for strongly correlated systems. These terms correspond to the summation of the infinite Born series.

T pp -G1-G2 equations for the Hubbard model
We now apply this result to the Hubbard Hamiltonian and find, where we introduced the particle-particle ladder term In the present case there exists only one distinct spin combination (two when considering ↑ ↔ ↓) of the particle pair that enters the single-particle EOM [cf. Eqs. (37) and (38)] which simplifies the equations. Numerical results for the T pp -G1-G2 scheme are presented in Sec. VII.
With that, the equation of motion for the time-diagonal two-particle Green function becomes, where the momentum representation of the particleparticle ladder term is given by For the T matrix in the particle-hole channel [28], the derivations of the single-time equations are performed in similar fashion as for the particle-particle T matrix in Sec. V A. The detailed derivation is given in Appendix B. Here, we summarize the main findings.

T ph -G1-G2 equations for a general basis
As for the GW and the TPP approximations, two quasi-Hamiltonians are introduced, and combined into a single quantity, The corresponding Møller operator of the particle-hole T matrix again obeys a time-dependent two-particle Schrdinger equation, The time derivative of G in TPH approximation follows as Again, for practical use, it is convenient to separate the correlation contributions from the mean-field terms via the introduction of an additional quantity: where the particle-hole ladder term is defined by As in the case of the particle-particle T matrix, Sec. V A, neglect of the Λ-terms exactly recovers the equation of motion for G in second-order Born approximation. Inclusion of theses terms, on the other hand, accounts for the entire Born series.

T ph -G1-G2 equations for the Hubbard basis
For the Hubbard system (for the definitions, see Sec. III D), we find, where we introduced the particle-hole ladder term for the Hubbard system Similar to the behavior in the TPP case, only one spin combination (two when considering ↑ ↔ ↓) contributes to the single-particle EOM in Eqs. (37) and (38). The T ph -G1-G2 scheme for the Hubbard model is numerically tested in Sec. VII.
With that, the equation of motion for the time-diagonal two-particle Green function becomes, with the momentum representation of the particle-hole ladder term, given by VI. DYNAMICALLY-SCREENED-LADDER APPROXIMATION So far we have considered three important selfenergy approximations: the second-Born approximation, GW and the particle-particle and particle-hole T matrices. While GW describes dynamical screening, for weakly coupled systems, the T -matrix selfenergy accounts for strong coupling but neglects dynamic screening effects. Therefore, the question arises how to combine strong coupling and dynamical screening into a single model in a computationally feasible way. An approximate to realize this within NEGF theory is the fluctuating-exchange approximation (FLEX) that combines T matrix and GW contributions according to Σ = Σ TPP + Σ TPH + Σ GW − 2Σ SOA , where the last term is needed to avoid double counting, for more details, see Ref. [28]. A fully selfconsistent treatment of dynamical-screening and strong-coupling effects is provided by the dynamically-screened-ladder approximation that has been studied in the context of the bound-state problem in a plasma medium in equilibrium [62]. For more details, see Ref. [63].
The G1-G2 scheme allows for a straightforward way to combine the GW (including exchange) and both Tmatrix approximations in a selfconsistent way for arbitrary nonequilibrium situations. This is achieved by including in the EOM of the time-diagonal two-particle Green function the terms with all effective Hamiltonians that were derived for GW , the particle-particle and the particle-hole T matrix, respectively, cf. Eqs. (56), (74) and (81). Then, the EOM for G, in a general basis becomes, Alternatively, we can rewrite this equation by using the polarization (Π) and ladder (Λ) terms that were defined by Eqs. (61), (77) and (84), where we combined both ladder terms into Obviously, Eq. (88) is a generalization of all previous cases: it additively includes the contributions of the second-order Born selfenergy (second line), polarization terms that account for dynamical screening and strong coupling terms. The SOA term that appears in each of the different approximations is included only once, so no double counting occurs. Since all contributions are treated on the same footing, this equation amounts to a simultaneous full account of dynamical screening and strong binary correlations. Alternatively, this approximation can be obtained from reduced-density-operator theory by neglecting threeparticle and higher correlations [5]; an early discussion was presented by Wang and Cassing [64]. It is easily verified that the entire Eq. (87) requires a CPU-time that has the same linear scaling with N t as all the special cases that were studied before. On the other hand, the polarization and ladder terms determine the scaling with the basis size N b . This is summarized in Tab. I and discussed in more detail in Sec. VII.

VII. VERIFICATION OF THE NUMERICAL SCALING
As was shown in the previous sections, the G1-G2 scheme reduces the time-diagonal Keldysh-Kadanoff-Baym equation within the HF-GKBA to a memory-less, time-local form. This means, the theoretical scaling is first order in the propagation duration. This dramatic acceleration is achieved by propagating, in addition to the single-particle Green function, also the time-diagonal two-particle Green function G. This function has, in general, four basis indices and, thus, a dimensionality of N 4 b , where N b is the single-particle basis dimension. The total scaling of the G1-G2 scheme with N b depends on the selfenergy and on the type of basis. In the following, we investigate this scaling more in detail, extending the analysis of Ref. [39].

A. Second-order Born selfenergy
We start by analyzing the N b -scaling of the SOA-HF-GKBA equation for G, Eq. (32), which we rewrite in a different form The r.h.s. of this equation contains four sums of dimensionality N b which are all independent of each other. They are evaluated by successive execution of the occurring tensor contractions. This means the total scaling of the CPU time, in this case, is of order N 5 b . For the Hubbard basis a first look at Eqs. (39) -(41) suggests an N 5 b -scaling, due to the commutator term in Eq. (39) and the summation in the Φ term of Eq. (41). However, in the Hubbard model the scaling can be further reduced. Note that the Hartree-Fock Hamiltonian, h HF (t), is a tridiagonal matrix and, thus, the commutator can be computed with N 4 b effort: On the other hand, the Φ term can be simplified by using the identity of Eq. (9): Here, the leading contribution to the difference, G < G < G < G < − G < G < G < G < , cancels (contribution with four functions G < ) which reduces the complexity. For the Hubbard basis, this reduces the numerical effort of the G1-G2 scheme to a N 4 b -scaling compared to the N 5 bscaling in the straightforward implementation [39]. In total, an acceleration is achieved for the SOA-G1-G2 scheme, compared to the ordinary HF-GKBA if N t N b , as summarized in Tab. I.
The reformulation above that eliminates products of four G < functions can be made for any basis choice. However, for the general basis this does not result in an improved N b -scaling. For the jellium basis the Eqs. (44) -(46) reveal a particularly favorable scaling with the basis size with N 3 b for which the above reformulation does not provide further improvement.

B. GW selfenergy
The additional terms of the GW approximation can change the N b scaling compared to the SOA case discussed in the previous section. For the general basis, the leadingorder terms for the scaling with the basis size are found in Eqs. (55) and (61)  For the Hubbard basis the polarization terms [Eqs. (64) and (65)] can be reformulated by again using Eq. (9) to get From this, it is obvious that, compared to the secondorder Born approximation, no further complexity is added for GW in the Hubbard case, and the scaling with the basis size remains N 4 b . To explore the N b -scaling for the jellium basis we recall the polarization term, Eq. (66), As one can see, the tensor contraction over k can be executed independently ofp. Thus, the full scaling of the GW -G1-G2 scheme for a jellium basis remains of order N 3 b , as in the case of the standard HF-GKBA.

C. T -matrix selfenergies
The T -matrix equations [Sec. V] behave very similar to the GW equations. For a general basis set with a four-index interaction tensor both, TPP and TPH scale as N 6 b which can be directly seen from Eqs. (73) and (77), as well as Eqs. (80) and (84).
For the Hubbard basis we can now use Eq. (9) to eliminate contributions that are of second order in G < from the ladder terms in Eq. (78), as well as in Eq. (85), For both cases one can see that the remaining scaling order of the equations is N 4 b since all internal summations have been eliminated.
In the jellium basis the T matrices show a different scaling behavior compared to GW . To see this, we reproduce the two ladder terms of Eqs. (79) and (86), Table I. Scaling of the CPU time with the number of time steps Nt and basis dimension N b of the traditional non-Markovian HF-GKBA and the present time-local scheme (G1-G2), for three relevant basis sets and the selfenergy approximations considered in this paper: the second-Born approximation (SOA), GW approximation (GW ), the particle-particle (TPP) and particle-hole (TPH) T matrices, and the dynamicallyscreened-ladder approximation (DSL). Last column: CPU speedup ratio of the G1-G2 scheme compared to standard HF-GKBA. For DSL, currently no standard HF-GKBA version exists. Note that full two-time NEGF simulations always have cubic scaling with Nt.
Evidently, in both cases the tensor contraction of k depends on all other momenta p,p, q. Thus, the final scaling with the basis size becomes of order N 4 b . A summary of the numerical scaling with the propagation duration and the basis size is presented in Tab. I. At the same time, any practical implementation of the G1-G2 scheme could, in principle, carry a large overhead that prevents to achieve the theoretical scaling with the simulation duration and the basis dimension within a relevant parameter range. We, therefore, have implemented the G1-G2 scheme for each of the selfenergies discussed in this paper and present representative numerical results in Sec. VII D.

D. Numerical results for the Hubbard basis
As we have shown above (cf. Tab. I), the Hubbard basis is the most unfavorable case for the G1-G2 scheme. Therefore, we choose this case for numerical demonstrations. In Ref. [39] we presented the first numerical tests of this scheme and demonstrated that, for finite Hubbard clusters the predicted linear scaling is indeed, achieved for SOA and GW selfenergies, already for rather small values N t . Here we extend these simulations to the T -matrix selfenergies and the DSL approximation. Furthermore, we explicitly verify the N b -scaling. As a first test, we verify that the derived formulas of the G1-G2 scheme are equivalent to the original (non-Markovian) HF-GKBA formulation. As a test case we consider, in Fig. 1 the time evolution in a Hubbard dimer for SOA, GW , TPP and TPH selfenergies. The agreement is excellent, and the deviations are mostly due to the original HF-GKBA, as discussed in Ref. [39].
Next, we verify the scaling with the basis dimension N b for the SOA selfenergy. In Fig. 2 we show results for a large number of Hubbard chains of varying length, N b = 2 . . . 100. We clearly confirm the N 5 b -scaling for the standard implementation of the G1-G2 scheme that uses Eq. (41) [39]. This asymptotic behavior is reached already for N b 20. The second curve is for the same setup but uses the optimization, Eq. (89). Again, the predicted improved scaling according to N 4 b is clearly identified, at least for N b 50. This confirms the expected speedup of the SOA-G1-G2 scheme compared to the standard HF-GKBA, if N t N b . Thus, even for the most unfavorable case of a Hubbard basis [cf. Tab. I] the scaling advantage should be reached already for small simulation durations. To explore the scaling with N t in more detail we have performed a series of simulations for all selfenergy approximations, comparing the standard HF-GKBA to the G1-G2 scheme. The results are shown in Fig. 3 and confirm the quadratic (cubic) scaling of the CPU time with N t , for the standard HF-GKBA with SOA (GW ) selfenergy. Similar cubic scaling is observed for the two T -matrix approximations (not shown) whereas simulations with DSL approximation are not possible, at the moment. Let us now turn to the G1-G2 results (dashed lines). Each of the curves exhibits the predicted linear scaling, already for N t 20. Interestingly, in the G1-G2 scheme, the CPU time required for the rather involved T -matrix approximations is only slightly above the time required for the comparatively simple SOA case. Equally remarkable is the observation that the GW and DSL approximations, which, in Hubbard, rely on cross-coupling spin components, are rather close to the former selfenergies.
Note that, for the present small system (10-site Hubbard chain) "break even" of the G1-G2 scheme is reached for all selfenergies compared to the ordinary SOA-HF-GKBA (dark blue curve) well below N t = 100 whereas the original GW -HF-GKBA (light blue) is unfavorable, practically from the start. For larger times, the ordinary GW -HF-GKBA quickly turns out unfeasible (e.g., for N t ∼ 10 3 it requires 10 4 times longer simulations than GW -G1-G2), and the same applies to the T -matrix selfenergies. Thus, we conclude that, it is not just a quantitative gain in CPU time that the G1-G2 scheme delivers but, in many cases, highly accurate simulations (beyond the simple SOA selfenergy) become possible at all that are (currently) impossible otherwise.  In particular, at increased coupling, U/J 2, SOA selfenergies are known to be inaccurate (for an analysis see Ref. [28]) and for reliable simulations, more advanced approximations are crucial. In that context the DSL approximation is particularly attractive because it contains the dominant correlation effects selfconsistently. Until now such simulations have only occasionally been reported, for very small systems and short propagation times. An example of a four-site Hubbard chain is shown in Fig. 4. We observe excellent agreement of our DSL-G1-G2 scheme to the Wang-Cassing approximation simulations of Akbari et al. [65] confirming the equivalence of the two approximations. The results show excellent quantitative agreement with exact diagonalization data (black curve), however, for times tJ/ 30 deviations are growing.

VIII. DISCUSSION AND OUTLOOK
In this paper we analyzed the properties of nonequilibrium Green functions in the frame of the generalized Kadanoff-Baym ansatz with Hartree-Fock propagators (HF-GKBA). Due to the non-Markovian structure of the collision integral, HF-GKBA simulations have an unfavorable quadratic (cubic) scaling with the number of time steps, for second-order Born (more complicated) selfenergies. At the same time, it has been reported earlier that this memory integral can be formally eliminated in favor of coupled time-local differential equations for the single-particle and two-particle density matrix [5,40]. An equivalent formulation in the framework of nonequilibrium Green functions has been established in Ref. [39]-the G1-G2 scheme. The formal equivalence between both approaches is important because it means that the G1-G2 scheme retains all attractive properties of the HF-GKBA: it is total-energy conserving and time-reversible [66]. Furthermore, all selfenergies from NEGF theory that have been derived, e.g. using diagrammatic techniques, can be transformed into a time-local form, by applying the HF-GKBA.
On the other hand, the former analyses concentrated, e.g., mainly on spatially homogeneous systems (jellium) [5] and did not include computational aspects such as the CPU time requirement. The scaling with the propagation time and basis size have only recently been analyzed in detail in conjunction with the G1-G2 scheme [39], and it was confirmed that the N 1 t -scaling can be achieved in practice. Here, we substantially extended these results, including additional high-level selfenergies such as the particle-particle and particle-hole T -matrix selfenergies and the screened-ladder approximation. In each case N 1 tscaling of the CPU time could be confirmed giving rise to a remarkable N 2 t -scaling advantage compared to the standard HF-GKBA scheme (Fig. 3) which was found to be independent of the single-particle basis used for the simulations. Furthermore, we re-analyzed the CPUtime scaling with the basis dimension N b and observed that the G1-G2 scheme has an overhead, compared to standard HF-GKBA, that is, at most, first order in N b , cf. Tab. I. Even for the most unfavorable basis-the Hubbard basis-the G1-G2 scheme has only a N 1 b overhead (down from a N 2 b overhead reported in Ref. [39]) which could be achieved by a reformulation of the scattering term in the G2-equation, cf. Sec. VII A. Thus, we expect that the G1-G2 scheme outperforms the standard HF-GKBA approach, in all cases of practical relevance, which can be seen from the CPU-time scaling ratio summarized in the right column of Tab. I.
With the G1-G2 scheme NEGF simulations (within the HF-GKBA) have been brought to the same CPU time scaling as many other time-dependent approaches, including semiclassical molecular dynamics, hydrodynamics, Boltzmann-type kinetic equations, TDDFT (adiabatic approximation), and the time-dependent Schrdinger equation. Most importantly, now long simulations are feasible that were previously prohibited by the memory structure (resulting in the N 2 t or N 3 t discussed above) without compromising the quality of the treatment of electronic correlations. We also showed that the inclusion of initial correlations in the G1-G2 scheme is trivial, and their propagation again requires a CPU time effort that is of order N t . Also the precomputation of the correlated initial state, e.g. via imaginary time stepping or adiabatic switching, see, e.g., Ref. [21], can be carried out separately and does not effect the propagation scaling.
While we presented numerical results only for the Hubbard model, even larger gains, compared to the standard HF-GKBA, are predicted for jellium (e.g. electron gas, dense quantum plasmas, electron-hole plasmas etc.) and for more general basis sets where the interaction tensor has four indices (e.g. electron dynamics in atoms and molecules). At the same time, the removal of the memory integral as the main CPU time bottleneck was achieved by computing the dynamics of an additional quantity-the time-diagonal two-particle Green function G ijkl . Thus, the new bottleneck in the G1-G2 scheme is the memory cost to store this four-dimensional tensor (only the current values are required), but this can be mitigated by suitable parallelization concepts.
By mapping NEGF simulations to a time-local scheme for single-time quantities, it should be expected that close connections exist with reduced-density-operator theory (RDO) [5,39,40]. The latter has been an independent many-body approach that has been successfully applied in many areas, including semiconductor optics, e.g. Refs. [67,68], dense plasmas [32], correlated electrons [54,65,69], nuclear matter [70], and cold atoms [71]. Our results indicate the correspondence between important selfenergy approximations of NEGF theory to closure relations of RDO and confirm and extend earlier results on the particle-particle T matrix [61] and the GW approximation [54]. We also investigated the simultaneous treatment of strong coupling and dynamical screening effects by combining ladder and polarization terms in the equation for G. This lead us to the dynamicallyscreened-ladder approximation (DSL), in Sec. VI. This approximation includes all two-particle interaction contributions and is, thus, equivalent to an approximation considered by Wang and Cassing before [64]. The equivalence of the two approximations was confirmed by the excellent agreement with the numerical results of Akbari et al. [65] for a small Hubbard cluster, cf. Fig. 4.
Despite the high quality of the DSL, we also observed that it is in quantitative agreement with exact diagonalization (CI) data (black curve in Fig. 4) only during the initial relaxation phase (for times tJ/ 30) [65]. So, clearly more systematic comparisons to CI results, for a broader range of coupling strengths and filling fractions, are desirable to understand the applicability limits of the DSL. While CI simulations are limited to very small particle numbers (basis size N b ) the G1-G2 scheme in DSL and simpler approximations can treat much larger systems. To go beyond those parameters where the DSL approximation is valid, further improved approximations are in high demand. This will require to partially include three-particle correlations. Examples are the Kirkwood superposition approximation of classical statistical physics [72] (for recent applications see Refs. [73,74]), the approximation by Nakatsuji and Yasuda [75,76], and selfenergy corrections to the BBGKY hierarchy [5]. Another route to improvements starts from nonequilibrium Green functions theory where one approach is to apply the GKBA but replace the Hartree-Fock propagators by correlated propagators [21]. Another concept is to replace the GKBA entirely, by an improved reconstruction ansatz. In both cases, the procedure outlined in the present paper will allow one to derive the corresponding improved G1-G2 scheme. Since the applicability limits of the GKBA are still not fully explored, full two-time NEGF simulations will remain indispensible for tests and benchmarks, see, e.g., Ref. [77].
In conclusion, let us come back to the remarkable capability of the G1-G2 scheme to efficiently perform long-time simulations of correlated electron dynamics. With this it should be feasible to reach thermodynamic equilibrium (or a quasi-stationary or pre-thermalized state) of the electrons. At the same time, slower processes, such as the equilibration with heavier particles (e.g. with the lattice in solids or with ions in dense plasmas) will make it desirable to develop a multiscale approach. This can be based on approximate solutions of the G1-G2 equations, e.g. by using retardation expansions [5] or the correlation-time approximation [78], eventually approaching the Markovian Boltzmann equation or local thermodynamic equilibrium. In that case a connection of the kinetic simulations to quantum hydrodynamic models, see, e.g., Refs. [79,80], could be a promising approach.

Symmetry relations
The single-particle time-evolution operator U fulfills the symmetry where G ij (t , t) has been used. Likewise, the two-particle propagator obeys, where Eq. (A1) has been used.

Group property
Utilizing Eqs. (22) and (23), now the group property for the propagator U is derived for all relevant time orderings. Starting with five different cases have to be considered. For t =t = t one gets k U ik (t, t)U kj (t, t) = k δ ik δ kj For t =t one gets k U ik (t, t)U kj (t, t ) = as well as fort = t , For t >t > t , the propagators reduce to U ij (t, t ) = G R ij (t, t ), for which Eq. (22) is directly applicable. For the analogous case, t <t < t , one obtains U ij (t, t ) = −G A ij (t, t ) which, together with Eq. (23), leads to, ijpq (t,t)U (2) pqkl (t, t ) , for the two-particle propagator. By construction, the particle-hole T matrix obeys the following symmetry [cf. Eq. (11)], T ph,≷ ijkl (t, t ) = T ph,≶ jilk (t , t) .