Numerical renormalization group method for computing four-point correlation functions

Four-point correlation functions commonly appear in various contexts of the theory of strongly correlated systems, including diagrammatic extensions of dynamical mean-field theory (DMFT). Here we develop the numerical renormalization group (NRG) method for computing four-point correlation functions in quantum impurity systems. First, we derive the Lehmann representation for general four-point functions (i) in imaginary Matsubara frequencies, (ii) on the real-frequency axes at zero temperature, and (iii) on the Keldysh contour. By using the complete basis of energy eigenstates constructed within NRG, four-point functions can be computed at arbitrarily low temperatures. We present results for paradigmatic models, including the effective quantum impurity model arising in DMFT treatments of the Hubbard model.