G1-G2 scheme: Dramatic acceleration of nonequilibrium Green functions simulations within the Hartree-Fock generalized Kadanoff-Baym ansatz

The time evolution in quantum many-body systems after external excitations is attracting high interest in many fields, including dense plasmas, correlated solids, laser-excited materials, or fermionic and bosonic atoms in optical lattices. 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) self energy, which has substantially extended the scope of NEGF simulations. In a recent Letter [Schlünzen et al., Phys. Rev. Lett. 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 the G1-G2 scheme. This allows one to perform the same simulations with order $T^1$ scaling, both for SOA and $GW$ self energies 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 self energies, including the $T$ matrix in the particle-particle and particle-hole channels and the dynamically-screened-ladder approximation. Finally, we demonstrate that, for all self energies, 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$, even for Hubbard systems.

Figure 1

Two examples of the Keldysh “round-trip” time contour that are used in NEGF theory to treat initial correlations, see, e.g., Ref. [2]. (a) Contour ${\mathcal{C}}_{A}S$ containing an initial real-time interval (from $-\infty$ to $t_0$) to adiabatically “switch on” the pair interaction, starting from an uncorrelated state. (b) Contour ${\mathcal{C}}_M$ with an imaginary branch allowing us to include thermodynamic-equilibrium correlations via Green functions which have one time argument on the real branch and one on the imaginary branch [33, 42].

Figure 2

Comparison of the ordinary HF-GKBA and the G1-G2 scheme for a Hubbard dimer with $U=J$ at half filling. The initial state was uncorrelated with the entire density concentrated at the first lattice site. Rows correspond to SOA, $GW$, TPP, and TPH self energies. Right column shows the deviation $\mathrm{\Delta }{n}_1\left(t\right)=n_1^{\mathrm{G}1-\mathrm{G}2}\left(t\right)-n_1^{\mathrm{ordinary}}\left(t\right)$ of the densities of both schemes on site 1.

Figure 3

CPU-time scaling of the SOA-G1-G2 scheme with the basis size $N_b$ comparing the direct, Eq. (39) [39], and the optimized implementation, Eq. (89). Results are for a 1D Hubbard chain.

Figure 4

CPU-time scaling with the simulation duration $N_t$, comparing the standard HF-GKBA to the G1-G2 scheme. G1-G2 data are shown for five self-energy approximations (indicated in the legend) all of which clearly exhibit linear scaling. In contrast, the standard HF-GKBA scales as $N_t^2$, for SOA, and $N_t^3$, for $GW$. Results are for a 10-site Hubbard chain.

Figure 5

Density evolution comparing the DSL-G1-G2 scheme to the results of Ref. [75] (WC) and exact-diagonalization simulations. The system is a four-site Hubbard chain at $U=0.1J$ and half filling; the simulations started from a noninteracting (uncorrelated) initial state, where the first two sites are doubly occupied.