Generalized Keldysh formalism for nonequilibrium correlation functions and its application to fluctuation dynamics

Abstract: Recent advances in time-resolved spectroscopies provide increasing access to collective dynamics in correlated quantum materials. However, computing the corresponding nonequilibrium two-particle correlation functions remains a major challenge. Here, by introducing a contour-dependent virtual probe field within the generalized Keldysh formalism, we propose an approach that computes such correlation functions with the vertex corrections essential for describing collective dynamics. In particular, we introduce a linear integral equation that computes the correlation functions without explicitly constructing the four-time vertex kernel, and develop its matrix-free Krylov solver based on quantics tensor trains. Combining our method with nonequilibrium dynamical mean-field theory, we show that the fluctuation dynamics of the order parameter in a nonequilibrium symmetry-broken state depends significantly on whether vertex corrections are included, and that the fluctuation and its decay time grow near the nonthermal critical point. Our approach thus provides a practical route for evaluating nonequilibrium correlation functions, which are emerging as key observables for characterizing states far from equilibrium.

• The paper introduces a generalized Keldysh approach to compute vertex-corrected nonequilibrium correlation functions using a virtual probe field.
• It leverages a quantics tensor train (QTT) matrix-free GMRES solver to reduce computational complexity while achieving high-resolution time evolution.
• Application to the Hubbard model shows that vertex corrections are essential for accurately capturing fluctuation dynamics and nonthermal criticality.

Generalized Keldysh Formalism for Nonequilibrium Correlation Functions and Fluctuation Dynamics

Motivation and Background

The computation of nonequilibrium two-particle correlation functions is a fundamental requirement for interpreting modern time-resolved spectroscopies, such as tr-Raman and time-resolved RIXS, which probe the real-time evolution of collective excitations in correlated quantum materials. The accurate evaluation of these correlation functions, including their vertex corrections, is crucial for understanding dynamical entanglement, relaxation, and criticality far from equilibrium. Traditional approaches based on nonequilibrium Green's function (NEGF) methods, as well as direct wave-function-based techniques, either fail to incorporate collective effects efficiently or become computationally prohibitive in large and finite-temperature systems. The Bethe-Salpeter equation (BSE) framework, while formally correct, is resource intensive due to the involvement of four-time kernels. Consequently, most practical calculations rely on bubble-level approximations, neglecting vertex corrections and, thereby, collective dynamics.

Formal Development: Generalized Keldysh Approach

This paper develops a unified formalism to compute vertex-corrected nonequilibrium correlation functions via the generalized Keldysh contour, leveraging a virtual probe field to extract the functional derivative response. The central advancement is the implementation of a contour-dependent Hamiltonian $\hat{H}_F = \hat{H} + F(z)\hat{A}$, where $F(z)$ is an external field introduced selectively on branches of the contour. This extension enables differentiation of the observable’s expectation value with respect to the probe field, yielding all physical components of the correlation function (Matsubara, left-mixing, retarded, lesser) beyond the equilibrium fluctuation-dissipation regime.

Figure 1: (a) The Keldysh contour $\mathcal{C}$ with forward, backward, and Matsubara branches, and the virtual probe field injection. (b) Diagrammatic representation of the linear integral equation for the response function. (c) Time-resolved fluctuation dynamics in a quench protocol.

A key theoretical insight is the formulation of a linear integral equation for the response function $G' = \delta G / \delta h$. Given a converged unperturbed Green's function $G$, the equation

$$G'(z_1, z_2; z_p) - [G * \Sigma'^{\rm int}[G, G'] * G](z_1, z_2; z_p) = \mathcal{B}(z_1, z_2; z_p)$$

is closed and linear in $G'$, with the kernel encoded in the action of the self-energy functional derivatives, circumventing the explicit construction of the four-time vertex kernel.

Numerical Implementation: Quantics Tensor Train and Matrix-Free Solver

The numerical solution employs a quantics tensor train (QTT) representation for the two-time functions $G$ and $G'$, enabling orders-of-magnitude compression and fine time resolution. The operator evaluation exploits tensor network contractions directly within the compressed format, allowing scalable computations (storage $\mathcal{O}(D_{\rm max}^3)$, operations $F(z)$0, where $F(z)$1 is the bond dimension). The linear equation is solved via GMRES with QTT-restarted cycles and dynamic reorthogonalization to stabilize basis truncation.

Algorithmic steps:

• Compute the probe-induced right-hand side via convolution.
• Apply the operator by contracting $F(z)$2, evaluating self-energy variations, and subtracting the sandwich term.
• Solve iteratively via restarted GMRES, maintaining a residual tolerance $F(z)$3.

Application to Nonequilibrium Fluctuation Dynamics

The methodology is demonstrated for order-parameter fluctuation dynamics in the half-filled Hubbard model on the Bethe lattice, following an interaction quench in antiferromagnetic initial conditions. The fluctuation $F(z)$4 traces the evolution of spin collective modes. Nonequilibrium DMFT is used as the embedding framework, with second-order perturbative impurity self-energies.

Figure 2: Quench dynamics of staggered magnetization $F(z)$5 and its fluctuation $F(z)$6 for different post-quench $F(z)$7 values; exponential vs. oscillatory behavior signals proximity to nonthermal critical points.

The results demonstrate that vertex corrections qualitatively alter the transient and long-time fluctuation behavior:

• With vertex corrections, the initial fluctuation response exhibits a decrease post-quench due to Hartree-Fock-driven Stoner suppression, followed by enhancement and slow decay.
• Without vertex corrections, the fluctuation growth is artificial and fails to capture collective relaxation effects.

  Figure 3: Critical dynamics: peak fluctuation, decay times, and Higgs-mode frequencies as functions of interaction strength, revealing the divergence and qualitative shift at nonthermal and thermal critical points.

Quantitative findings:

• Near the nonthermal critical point ($F(z)$8), both the maximum fluctuation amplitude $F(z)$9 and the decay time $\mathcal{C}$0 exhibit critical enhancement, indicating that two-particle observables, not just order parameters, encode nonequilibrium criticality.
• The difference between calculations with and without vertex corrections increases as the system approaches the critical regime.

Comparison and Benchmarking

Direct pump-probe simulations (subtracting Green's functions with and without probe field) and the linear-response perturbative approach yield consistent correlation functions at short times where numerical convergence is tractable. For longer times and broad time domains, the perturbative method retains stability and accuracy due to the matrix-free QTT-GMRES solver.

Figure 4: Comparison of spin-spin correlation functions computed via GCH pump-probe (solid lines) and perturbative GCH (points) in equilibrium and nonequilibrium, confirming methodological consistency.

Implications and Future Directions

Practically, the generalized Keldysh formalism and QTT-based linear solvers facilitate computation of full vertex-corrected correlation functions for large, high-dimensional systems, including the nine-component objects required for time-resolved spectra. The ability to characterize nonequilibrium entanglement, relaxation, and criticality from two-particle correlations enables direct connection to observables in pump-probe and quench experiments. Theoretical implications include rigorous access to quantum Fisher information and multipartite entanglement dynamics. Extensions to multi-orbital systems, realistic lattice geometries, and QMC-evaluated self-energies are feasible as the framework only relies on two-time operator action, permitting further scalability.

Figure 5: Equilibrium order parameter and fluctuation as functions of $\mathcal{C}$1, with the fluctuation diverging at the critical point.

Figure 6: Contour structure and probe field assignment, relating the probe protocol to physical correlation components.

Figure 7: Diagrammatic comparison between standard BSE and the perturbative GCH linear-response scheme.

Figure 8: Fit of Higgs-mode oscillation frequency $\mathcal{C}$2 from time-series analysis of $\mathcal{C}$3 across interaction regimes.

Conclusion

This work presents a robust formalism for the calculation of nonequilibrium vertex-corrected correlation functions within the generalized Keldysh approach, introducing a QTT-based matrix-free linear solver for efficient high-resolution computation. Application to fluctuation dynamics in the Hubbard model demonstrates that vertex corrections are indispensable for capturing collective relaxation and criticality, expanding the toolkit for theoretical interpretation of quantum materials in nonequilibrium regimes. The formulated approach is poised to impact future studies of multi-component entanglement, pump-probe spectroscopy, and emergent phenomena in correlated electron systems (2607.11055).