Nonequilibrium Green Functions Simulations for Large Correlated Systems

Abstract: Correlated real-time dynamics in large, spatially inhomogeneous quantum systems remain difficult to access with nonequilibrium many-body methods. Two-time nonequilibrium Green functions (NEGF) retain dynamical correlations but their computational runtime grows cubically with the number of time steps $N_\mathrm{t}$. This scaling bottleneck could recently be overcome by introducing the G1--G2 scheme that is linear in $N_\mathrm{t}$, but requires propagation of a two-particle correlation function and may suffer from numerical instabilities. This has restricted simulations to small systems with $N_\mathrm{b} \sim 10<sup>2$ basis states. Here we introduce a quantum-fluctuation formulation of nonequilibrium Green functions, denoted $δ$NEGF, that represents dynamical two-particle correlations through fluctuations of field-operator products, $δ\hat G$. This guarantees stable dynamics by preserving the positivity of the reduced density matrices, avoids the explicit storage of the two-particle Green function, and reduces the propagation to a finite ensemble of Hartree-Fock-like trajectories. Combined with a stochastic low-rank decomposition of the correlation functions, the method retains time-linear scaling while extending dynamical $GW$ and particle-particle and particle-hole $T$-matrix simulations to basis sizes of order $N_\mathrm{b}\sim 10<sup>4$. We benchmark $δ$NEGF against exact and HF-GKBA results for lattice systems, finding stable correlated dynamics also at strong coupling. We further demonstrate large-scale simulations of diffusion in two-dimensional Hubbard lattices and ultrafast relaxation in graphene nanoribbon heterostructures with long-range Coulomb interactions. These results establish $δ$NEGF as a scalable route to dynamical self-energy simulations of large, spatially inhomogeneous correlated quantum systems beyond the reach of existing NEGF implementations.

• The paper introduces δNEGF, which represents two-particle correlations as positive-semidefinite field fluctuations and reduces propagation to independently parallelized, low-rank trajectories.
• The method changes general runtime from cubic time-step scaling to approximately NₛNᵦ⁴Nₜ and supports simulations with basis sizes near 10⁴, while preserving key conservation properties in specified interaction models.
• Benchmarks show stable and accurate results for strong-coupling charge dynamics, large-lattice diffusion, spin correlations, Mott-insulator spectra, and nonequilibrium responses in graphene nanostructures.

Motivation and scope

Correlated real-time dynamics in large, spatially inhomogeneous quantum systems remain difficult to access with nonequilibrium many-body methods. Conventional two-time nonequilibrium Green functions (NEGF) retain dynamical correlations and memory, but their runtime scales cubically with the number of time steps $N_\mathrm{t}$, and correlated calculations have historically been restricted to small model systems or symmetry-reduced settings. The G1–G2 scheme—a time-local reformulation of the HF-GKBA—removed the cubic time scaling, but it requires storage and propagation of the two-particle correlation function $\mathcal{G}$, a rank-4 tensor with $N_\mathrm{b}^4$ elements, and inherits numerical instabilities from Hartree–Fock propagators. In practice, this has confined dynamical NEGF simulations to basis sizes of order $N_\mathrm{b}\sim 10^2$.

This paper introduces a quantum-fluctuation formulation of NEGF theory, denoted $\delta$NEGF, that removes both bottlenecks simultaneously: two-particle correlations are represented through fluctuations of field-operator products, $\delta\hat{G}$, whose factorized form guarantees positive semi-definiteness of the reduced density matrices and permits a stochastic low-rank decomposition. The result is an ensemble of Hartree–Fock-like fluctuation trajectories that can be propagated independently and in parallel, extending dynamical $GW$, particle-particle ($\delta$TPP), and particle-hole ($\delta$TPH) $T$-matrix simulations to basis sizes of order $\mathcal{G}$0.

Theoretical framework

The method builds on the channel decomposition of self-energy approximations into particle-particle (pp), longitudinal particle-hole (ph), and transversal particle-hole ($\mathcal{G}$1) channels, associated respectively with the TPP, TPH(+X), and $\mathcal{G}$2(+X) approximations. The central object is the generalized susceptibility $\mathcal{G}$3, which satisfies a Bethe-Salpeter equation on the Keldysh contour. The key observation is that each susceptibility can be written as a contour-ordered correlation function of single-particle fluctuation operators,

$\mathcal{G}$4

where $\mathcal{G}$5 is defined as the deviation of a two-operator product from its expectation value. Solving equations of motion for these fluctuations—rather than for two-particle quantities directly—yields several structural advantages.

Positivity by construction: because $\mathcal{G}$6 is expressed as expectation values of fluctuation operators and their adjoints, it is positive semi-definite at all times. This preserves the $\mathcal{G}$7-representability of the corresponding reduced density matrices and eliminates the instabilities that plague HF-GKBA propagations at strong coupling, without purification or other regularization procedures.

Linearization hierarchy: three levels of approximation are introduced for the particle-hole channels—the approximation of second moments (yielding $\mathcal{G}$8RPA), the polarization approximation (PA), and the reduced polarization approximation (yielding $\mathcal{G}$9 and $N_\mathrm{b}^4$0TPH)—plus a ladder approximation (LA) for the pp channel (yielding $N_\mathrm{b}^4$1TPP). The authors show that these correspond to an intermediate level of self-consistency between the HF-GKBA and a full mean-field treatment of the Bethe-Salpeter equation.

Conservation laws: the Keldysh component satisfies the exchange symmetry $N_\mathrm{b}^4$2 in all channels, which is the necessary condition for conserving behavior. The paper demonstrates analytically and numerically that all approximations used are energy conserving for local Hubbard-type interactions, with $N_\mathrm{b}^4$3RPA and $N_\mathrm{b}^4$4+X additionally conserving for diagonal (Coulomb-type) interactions.

Scaling and stochastic low-rank decomposition

The decisive computational advantage comes from combining the fluctuation representation with a low-rank factorization of the initial state. Since only the initial susceptibility needs to be decomposed—and the linearized dynamics preserve rank throughout the evolution—the problem reduces to propagating $N_\mathrm{b}^4$5 independent single-particle fluctuation trajectories. The scaling comparison is stark:

Basis	Standard KBE	G1–G2	$N_\mathrm{b}^4$6NEGF
Runtime (general)	$N_\mathrm{b}^4$7	$N_\mathrm{b}^4$8	$N_\mathrm{b}^4$9
Runtime (Hubbard)	$N_\mathrm{b}\sim 10^2$0	$N_\mathrm{b}\sim 10^2$1	$N_\mathrm{b}\sim 10^2$2
Memory (general)	$N_\mathrm{b}\sim 10^2$3	$N_\mathrm{b}\sim 10^2$4	$N_\mathrm{b}\sim 10^2$5
Memory (Hubbard)	$N_\mathrm{b}\sim 10^2$6	$N_\mathrm{b}\sim 10^2$7	$N_\mathrm{b}\sim 10^2$8

The initial state can be generated efficiently by drawing random samples from a Gaussian distribution with covariance given by the ideal susceptibility; by the central limit theorem, the reconstruction error scales as $N_\mathrm{b}\sim 10^2$9 and is effectively independent of system size. A further distinguishing feature relative to HODLR- or QTT-based compression schemes is that the difference between approximate and exact rank remains constant in time, since the employed approximations are rank-preserving by construction.

A notable practical benefit is direct access to two-time response and correlation functions—density, spin, current, dipole, and pair fluctuations—from a single-time propagation, enabling calculation of optical conductivity, dynamic structure factors, and screened interactions without explicit two-time-plane propagation.

Numerical benchmarks

The paper validates the approach across five systems spanning weak to strong coupling, using DMRG, exact diagonalization, and dynamical DMRG as references.

CDW melting in a 1D Hubbard chain ($\delta$0TPP): For an 11-site charge-density-wave system at strong coupling $\delta$1, GKBA+TPP becomes unstable and breaks down for $\delta$2, while GKBA+2B fails entirely. In contrast, $\delta$3TPP remains stable over the entire coupling range $\delta$4 to $\delta$5 and reproduces the exact density imbalance and correlated double occupation with high accuracy. This constitutes the strongest claim of the paper: stable, quantitatively accurate correlated dynamics at coupling strengths previously inaccessible to any NEGF-based simulation.

Diffusion in 2D Hubbard lattices ($\delta$6TPP): For a $\delta$7 lattice with 74 particles, convergence studies show that global observables such as the expansion velocity are accurately captured with as few as 32 samples—a speedup exceeding $\delta$8—while space-resolved densities require roughly $\delta$9 samples. Scaling up to a $\delta\hat{G}$0 lattice with 2514 particles (using $\delta\hat{G}$1 instead of the ~81.6 million samples required for exact decomposition), the simulations reveal a qualitatively new phenomenon absent in smaller systems: at moderate coupling $\delta\hat{G}$2, expansion proceeds preferentially along the lattice diagonals, with a velocity ratio approaching 4, emerging only when correlations (TPP-level), sufficiently strong coupling, and sufficiently large system size ($\delta\hat{G}$3) are all present simultaneously. This demonstrates that finite-size effects in previous simulations masked genuine large-system physics.

Spin correlations in a Hubbard ladder ($\delta\hat{G}$4TPH): For a $\delta\hat{G}$5 ladder, the 2B approximation fundamentally fails to capture nonlocal antiferromagnetic correlations, and GKBA+TPH unphysically suppresses leg correlations beyond $\delta\hat{G}$6. The $\delta\hat{G}$7TPH approximation correctly captures the qualitative growth of both leg and rung correlations in the strongly correlated regime, though it quantitatively overestimates their magnitudes. In the nonequilibrium laser-driven case, both TPH variants capture the lower-energy spectral features missed by 2B but fail to reproduce the highest-energy peak near $\delta\hat{G}$8.

Optical response of a Mott insulator ($\delta\hat{G}$9): For a 128-site half-filled Hubbard chain, $GW$0 reproduces the qualitative opening of the Mott gap and correctly captures secondary peak structures, though it overestimates the gap and underestimates spectral broadening compared to DDMRG reference data. In the time-resolved dynamic structure factor following a laser pulse on a 12-site chain, transient Floquet-like sidebands and a persistent low-energy charge-excitation band are captured well, but the gapped Mott branch is over-bleached relative to exact diagonalization.

Carbon nanostructures with long-range Coulomb interactions ($GW$1RPA, $GW$2+X): For benzene and naphthalene within the extended Hubbard (Pariser-Parr-Pople) model, equilibrium spectra coincide with TDH/RPA and TDHF/RPAx respectively—but critically, this agreement breaks down under nonequilibrium driving, where TDH/TDHF strongly underestimate laser-induced spectral redistribution while the self-consistent $GW$3NEGF treatment captures bleaching and photoinduced structures. An application to a 7–9 armchair graphene nanoribbon heterostructure with 1728 sites resolves in-gap edge and junction states and reveals a delayed (~10 fs) build-up of excitonic population after pulsed excitation, consistent with earlier results on smaller systems.

Limitations and open questions

The paper is candid about several restrictions. First, the $GW$4NEGF approximations operate at a reduced level of self-consistency relative to the HF-GKBA, which can lead to a delayed onset of short-time relaxation dynamics—a systematic feature attributed to the neglected $GW$5 term in the susceptibility equation of motion. Second, energy conservation holds only for specific interaction classes: $GW$6TPP is conserving for local (Hubbard-type) interactions, while $GW$7RPA and $GW$8+X extend to diagonal (Coulomb-type) interactions, but not generally. Third, the accuracy of two-particle spectra at the $GW$9 level is comparable to RPA-like results, with known deficiencies such as gap overestimation and over-bleaching of gapped excitations. Fourth, the favorable scaling assumes factorizable interaction tensors; plane-wave bases for uniform systems remain less favorable and require additional symmetry exploitation. Finally, correlated single-particle spectral functions cannot be directly extracted within the current framework, and extension to open quantum systems, dynamically screened ladder self-energies, and embedding schemes for transport remains open.

Conclusion

This work establishes $\delta$0NEGF as a scalable formulation of dynamical self-energy simulations that simultaneously achieves time-linear propagation, guaranteed positivity of reduced density matrices, and basis sizes two orders of magnitude beyond previous dynamical NEGF calculations. The combination of rigorous benchmarking against exact methods—from CDW melting at $\delta$1 to spin correlations and Mott-gap spectroscopy—with demonstrably converged large-scale applications to $\delta$2-site Hubbard lattices and graphene nanoribbon heterostructures provides a credible foundation for treating heterogeneous correlated systems, including defects, interfaces, and molecular junctions, that were previously accessible only to adiabatic or mean-field descriptions.