Multipoint correlation functions: spectral representation and numerical evaluation

The many-body problem is usually approached from one of two perspectives: the first originates from an action and is based on Feynman diagrams, the second is centered around a Hamiltonian and deals with quantum states and operators. The connection between results obtained in either way is made through spectral (or Lehmann) representations, well known for two-point correlation functions. Here, we complete this picture by deriving generalized spectral representations for multipoint correlation functions that apply in all of the commonly used many-body frameworks: the imaginary-frequency Matsubara and the real-frequency zero-temperature and Keldysh formalisms. Our approach separates spectral from time-ordering properties and thereby elucidates the relation between the three formalisms. The spectral representations of multipoint correlation functions consist of partial spectral functions and convolution kernels. The former are formalism independent but system specific; the latter are system independent but formalism specific. Using a numerical renormalization group (NRG) method described in the accompanying paper, we present numerical results for selected quantum impurity models. We focus on the four-point vertex (effective interaction) obtained for the single-impurity Anderson model and for the dynamical mean-field theory (DMFT) solution of the one-band Hubbard model. In the Matsubara formalism, we analyze the evolution of the vertex down to very low temperatures and describe the crossover from strongly interacting particles to weakly interacting quasiparticles. In the Keldysh formalism, we first benchmark our results at weak and infinitely strong interaction and then reveal the rich real-frequency structure of the DMFT vertex in the coexistence regime of a metallic and insulating solution.


I. INTRODUCTION
A. Multipoint correlation functions A major element in the ongoing challenge of the quantum many-body problem is to extend our understanding, our analytical and numerical control, from the single-to the many-particle level. One-particle correlation functions, describing the propagation of a particle in a potentially highly complex, interacting environment, are clearly important. However, almost since the beginning of interest in the many-body problem, two-particle correlation functions have played an equally important role. They describe the effective interaction between two particles in the many-body environment, response functions to optical or magnetic probes, collective modes, bound states, and pairing instabilities, to name but a few, and are essential ingredients in Landau's Fermi-liquid theory [1].
For a long time, two-particle or four-point (4p) functions could only be computed by perturbative means, later including resummation and renormalization group schemes as well. Recently, the advance of numerical techniques to compute such functions fully nonperturbatively, albeit only locally, has opened a new chapter. For various impurity models (general ones as well as those arising in dynamical mean-field theory (DMFT) [2]), it was even found that 4p functions can exhibit divergences in * These authors contributed equally to this work. strongly correlated regimes [3][4][5][6][7][8][9][10][11][12][13]. These divergences occur for a special class of 4p functions, namely two-particle irreducible vertices, and have renewed the interest in thoroughly understanding the properties of 4p functions.
The properties of multi-or -point ( p) functions depend on the theoretical framework employed. Technically, p functions are correlation functions of operators with arguments. If the operators are taken at different times, the Fourier-transformed function depends on frequencies.
Clearly, these times and frequencies are real numbers, and the early works directly dealt with these real frequencies.
Later, the many-body theory was revolutionized by the Matsubara technique, based on a Wick rotation from real to imaginary times [14]. This significantly facilitates numerical approaches, and the vast majority of numerical work on p functions is based on the Matsubara formalism (MF). The periodicity in imaginary times yields functions of discrete imaginary or Matsubara frequencies. Physical results are obtained in the end by performing an analytic continuation, back from imaginary to real frequencies. Yet, crucially, the numerical analytic continuation from imaginary to real frequencies is an ill-conditioned problem [15]; though it may work fairly well for docile 2p functions, it is unfeasible for higher-point objects.
Formally, real-frequency frameworks are as well established as the MF: the zero-temperature formalism (ZF) works with ground-state expectation values of timeordered operators [14]. Finite-temperature real-frequency dynamics, both in and out of equilibrium, are treated in the Keldysh formalism (KF), based on a doubled time contour [16,17]. Due to the problematic analytic continuation and growing interest in nonequilibrium dynamics, the KF has gained much popularity in recent years. However, there are hardly any numerical real-frequency results of p functions depending on more than one frequency. Conceptually, it is therefore of great interest to extend the understanding of p functions that has been reached for imaginary frequencies to the real(-frequency) world.

B. Our approach
In this work, we achieve a deepened understanding of p functions, applying naturally in any of the real-or imaginary-time frameworks. The central principle of our approach is to separate properties due to time ordering from spectral properties of the system. The former govern the analytic structure of p functions; the latter are the key objects of the numerical evaluation. While exploiting this separation is standard practice for 2p functions, to our knowledge it has not yet been pursued systematically for p functions. Here, we derive spectral representations of general p functions for arbitrary , in all three frameworks (ZF, MF, KF). We illustrate the power of our approach by presenting numerical results for various local 4p examples.
The spectral representations are based on a density matrix, which defines the quantum averages, and an expansion in eigenstates of the Hamiltonian, which determines the time evolution. The spectral representations have similar forms in all three frameworks. Indeed, the spectral or Lehmann representation of 2p functions is arguably the most transparent way of demonstrating the analytic relation between Matsubara and retarded propagators [18]. Here, we achieve a similar level of transparency: all p functions are expressed by a sum over formalismindependent "partial" spectral functions, convolved with formalism-specific kernels. The partial spectral functions (PSFs) serve as unique and compact porters of the systemspecific information. Convolving them with the rather simple and system-independent kernels yields correlation functions in the desired form. Through this convolution, the correlation functions acquire their shape and characteristic long tails. In fact, the numerical storage of p functions, depending on −1 frequencies, can be problematic [19][20][21][22][23][24]. In this regard, the PSFs are advantageous since they have compact support, bounded by the largest energy scale in the system. For a set of external time arguments of a correlation function G(t 1 ,..., t ), the time-ordering prescription determines the corresponding order of operators in the expectation value. In our approach, we first consider all possible operator orderings, i.e., O 1 (t 1 ) · · · O (t ) for all permutations p(1, ... , ) = (1, ... , ). The Fourier transform of each of these ! expectation values yields a PSF; the collection of all PSFs describes the spectrum of the system (on the p level). Second, for given external times, some operator orderings are specified, depending on the formalism. In the ZF and MF, this is always one ordering; in the KF using the Keldysh basis, it may be a linear combination of multiple ones. The ordering is specified by kernel functions which are multiplicative in the time domain. In the frequency domain, they become convolution kernels, mapping PSFs to contributions to G(ω 1 ,..., ω ).
The spectral representations clearly reveal the similarities and differences of the various p functions. In the most complex setup of the KF with Keldysh basis, where each argument has an extra Keldysh index 1 or 2, we find that those components with a single Keldysh index equal 2 at position η, dubbed G [η] , have the simplest structure. Indeed, they are the (fully) retarded objects which can be obtained from Matsubara p functions via a suitable analytic continuation, iω i → ω i ± i0 + [25]. This does not apply to the remaining Keldysh components. However, we find that their convolution kernels K are linear combinations of the retarded kernels K [η] . Thereby, our spectral representations offer a direct way of discussing the relation between the Matsubara and all Keldysh p functions in explicit detail-this will be the topic of a forthcoming publication [26].
To illustrate our approach, we numerically evaluate the spectral representation of local 4p functions for selected quantum impurity models via the numerical renormalization group (NRG) [27]. Since its invention by Wilson in 1975 [28], NRG has become the gold standard for solving impurity models. It has the unique advantage of allowing one to (i) directly compute real-and imaginary-frequency results without the need for analytic continuation, (ii) reach arbitrarily low temperatures with marginal increase in the numerical costs, (iii) access vastly different energy scales, zooming in on the lowest excitation energies. Our NRG scheme is based on the full density matrix (fdm) NRG [29,30] in an efficient tensor-network formulation [31][32][33], with an additional, iterative structure to finely resolve regimes of frequencies |ω i | |ω j |, i = j. The rather intricate prescription for how to do this for 3p and 4p functions is explained in the accompanying paper [34].

C. Structure of this paper
Our spectral representations of p functions and the NRG scheme to evaluate them have numerous potential applications, such as finite-temperature formulations of Fermi-liquid theory and nonlocal extensions of DMFT in the MF at low temperatures or in real-frequency KF implementations. We elaborate on the applications in the concluding Sec. V B.
The rest of the paper is organized as follows. In Sec. II, we derive the spectral representations in a general setting, for arbitrary , a given density matrix, and any timeindependent Hamiltonian. After motivating our approach and introducing notation for = 2, we subsequently derive our results in the ZF, MF, and KF, the latter both in the contour and Keldysh bases. The main results of this section are the spectral representations, involving the PSFs (28) and summarized in Eqs. (26) and (27) for the ZF, in Eqs. (39) and (45) or (46) for the MF, and in Eqs. (63), (64) or (67) for the KF in the Keldysh basis.
In Sec. III, we briefly describe the quantities of interest for the numerical evaluation, i.e., local 4p correlation and vertex functions. Our numerical results are contained in Sec. IV. We start with a simple model for x-ray absorption in metals treated in the ZF and analyze power laws in the 4p vertex. Proceeding with the Anderson impurity model (AIM) in the MF, we first benchmark our results at intermediate temperatures against Monte Carlo data and then extend these results to lower temperatures to enter the Fermi-liquid regime and deduce the quasiparticle interaction. Thereafter, we treat the AIM in the KF, first testing our method at weak and infinitely strong interaction, before moving to the intermediate, strongly interacting regime. Finally, we present results for the DMFT solution of the one-band Hubbard model and compare the MF and KF vertex for both the metallic and insulating solution. In Sec. V, we summarize our results and give an outlook on applications. Appendices A to C are devoted to exemplary calculations in each of the ZF, MF, and KF. Appendix D describes details needed for amputating external legs when computing the ZF or KF 4p vertex, and App. E discusses an example of anomalous parts in both the MF and KF.

A. Motivation of partial spectral functions
To set the stage and introduce notation, we review the standard derivation of spectral representations for 2p correlators. We denote complete sets of energy eigenstates by underlined integers, e.g. {|1 }, with eigenvalues E 1 . We use calligraphic or roman symbols for operators or their matrix elements, A 12 = 1|A|2 . Expectation values are obtained through the density matrix , which, in thermal equilibrium at temperature T = 1/β, directly follows from the Hamiltonian H: For instance, we have AB = 12 ρ 1 A 12 B 21 , with ρ 1 = 1| |1 = e −βE1 /Z. The ZF assumes ρ 1 = δ 1g at T = 0 with a nondegenerate ground state |g [35]; the MF and KF work at any T . Further, the KF can also be used with a nonequilibrium density matrix. Yet, for simplicity, we do not consider such cases explicitly. In the ZF or MF, operators obey Hamiltonian evolution in real time t or imaginary time τ , respectively: The corresponding time-ordered (T ) correlators and their Fourier transforms in the ZF or MF are defined as where, depending on context, ω denotes a continuous real frequency or a discrete Matsubara frequency. (For brevity, we distinguish G in the ZF and MF solely through its arguments, t vs. τ or ω vs. iω.) Textbook calculations, based on judicious insertions of the identity in the form 1 = 1 |1 1| = 2 |2 2|, yield the following Lehmann representation of the Fourier-transformed ZF correlator: Here, we used ζ = 1 (ζ = −1) for bosonic (fermionic) operators, ω ± = ω ± i0 + for convergence of real-time integrals, and E 21 = E 2 −E 1 . The MF correlator is obtained as In the bosonic case, where iω can equal zero, terms with both iω = 0 and E 21 = 0 (if present) can be (analytically) dealt with by taking the limit E 21 → 0 in Eq. (5b), yielding (1 − e −βE21 )/E 21 → β. It will sometimes be convenient to make this "anomalous" case explicit, writing We can also consider 2p correlators arising in the KF, e.g. G +− (t) = −i A(t)B , or the retarded propagator G 21 (t) = −iθ(t) [A(t), B] −ζ , where θ denotes the step function, [·, ·] − a commutator, and [·, ·] + an anticommutator. Their Fourier transforms read Evidently, the imaginary-time and the retarded correlator are connected by the well-known analytic continuation, G 21 (ω) = G(iω → ω + ) [36]. They are often expressed through the "standard" spectral function S std : It would be convenient to have spectral representations capable of describing G(ω) and G +− (ω), too. To this end, we define the partial spectral functions (PSFs) Clearly, we can reconstruct S std from S according to Furthermore, S can be used to express all correlators encountered so far: For the bosonic case, the representations (8b) and (11b) involve a subtlety: To correctly reproduce the anomalous term of Eq. (6), the integral over ω must be performed first and the limit E 21 → 0 taken only thereafter. If, instead, S std (ω) is simplified first by using E 21 = 0 to conclude that ρ 1 − ρ 2 = 0, the anomalous terms are missed. To ensure that Eq. (6) is always correctly reproduced, including its anomalous terms, we refine the representation (11b) by using a kernel in which the anomalous case iω − ω = 0 is specified separately: The sum in Eq. (12a) is over the two permutations of two indices, p = (1 2) ∈ {(1 2), (2 1)}, with corresponding signs ζ (12) = 1, ζ (21) = ζ. Furthermore, O 1 , ω 1 , and ω 1 are components of operator and frequency tuples defined and (ω 1 , ω 2 ), respectively, the latter being (dummy) integration variables. Spelled out, the representation (12) yields which matches Eq. (6) after relabeling 1 ↔ 2 in the second term. The anomalous part, relevant only for ζ = 1 and ω = 0, receives equal contributions from both terms, since E 21 = 0 implies ρ 1 = ρ 2 . For p functions, the treatment of anomalous terms is more involved, as bosonic frequencies can arise through several combinations of fermionic ones. The formalism to be developed below will enable us to deal with such subtleties in a systematic fashion. In the following sections, we derive spectral representations for p functions. We start with the ZF, where the derivation is most compact. Subsequently, we show that the same PSFs also arise in the MF. Finally, we extend our analysis to the KF in the contour and Keldysh bases.

B. Zero-temperature formalism
To obtain concise expressions for p functions, we introduce further compact notation. A priori, such functions of operators O = (O 1 ,..., O ) depend on time or frequency arguments, t = (t 1 ,..., t ) or ω = (ω 1 ,..., ω ). Permutations of such ordered tuples are denoted O p = (O 1 ,..., O ), t p = (t 1 ,..., t ), and ω p = (ω 1 ,..., ω ). Here, by definition, the permutation p(12... ) = (12... ) [or p = (12... ) for short] acts on the index tuple (12... ) by replacing i by p(i) = i in slot i. Note that p moves i to slot j = p −1 (i), replacing j there by p(j) = i. If p = (312), e.g., then (t 1 , t 2 , t 3 ) p = (t 3 , t 1 , t 2 ). We also use the shorthands Because of time-translation invariance, p functions actually depend on only − 1 independent time or frequency arguments. Nevertheless, as anticipated above, it will be helpful to use all frequencies, it being understood that their sum equals zero. With the shorthand ω i···j = j n=i ω n , we can express the corresponding energyconservation relations as We use calligraphic symbols, G, K, S, for functions of all arguments and roman symbols, G, K, S, for functions of the independent −1 arguments. In the time domain, we define the time-ordered ZF correlator as It is invariant under a uniform shift of all times, e.g. t i → t i −t , and thus depends on only −1 time differences, (t 1 − t ,..., t −1 − t , 0). Accordingly, in the frequency domain, We write G(ω), the part remaining after factoring out 2πδ(ω 1··· ), with frequency arguments, it being understood that they satisfy energy conservation, ω 1··· = 0. The effect of the time-ordering procedure in the definition of G(t) can be expressed as a sum over permutations involving products of step functions θ: For a given -tuple of times t = (t1,..., t ), only that permutation p in Eq. (17) yields a nonzero contribution, K(tp) = 0, for which the permuted times, tp = (t 1 ,..., t ), satisfy t 1 > t 2 > · · · > t . In the corresponding operator product, S[Op](tp), larger times appear to the left of smaller ones. Thus, the right-to-left order of operators in the product matches the order in which their times appear on the time axis, drawn with times increasing toward the left.
In Eq. (17a), the sum is over all permutations p of the indices labeling operators and times. For a given choice of times, only one permutation yields a nonzero result, namely the one which arranges the operators in timeordered fashion (larger times left of smaller times, cf. Fig. 1). The sign ζ p is +1 (−1) if O p differs from O by an even (odd) number of transpositions of fermionic operators. Our Eqs. (17) conveniently separate properties due to time ordering, contained in the kernel K, from spectral properties involving eigenenergies and matrix elements, contained in S. This separation is well known for = 2; for larger , it was pointed out in 1962 by Kobe [37], but without elaborating its consequences in great detail. The guiding principle of this paper is to systematically exploit this separation-on the one hand, to unravel the analytic structure of p correlators, arising from K, on the other hand, to facilitate their numerical computation, involving mainly S. For notational simplicity, we start by considering the identity permutation, p = id, in Eq. (17); generalizing to arbitrary p afterward will be straightforward. The multiplicative structure of Eq. (17a) translates into an -fold convolution in frequency space: We may thus Fourier transform K and S separately. We begin with the Fourier transform of K: It is convenient to use t i = t i − t i+1 for i < and t = t as independent integration variables. Then, t i = t i··· (short for j=i t j ) and ω · t = i=1 ω 1···i t i , such that The frequencies ω + 1···i = ω 1···i + i0 + again contain infinitesimal imaginary parts to ensure convergence of real-time integrals. As for G(ω), the arguments of K(ω) are understood to satisfy energy conservation, ω 1··· = 0. A similar result, with imaginary parts sent to zero, was reached by Kobe [37]. Importantly, however, for numerical calculations, the imaginary parts are necessarily finite and their mutual relations must be treated carefully-this is discussed in Sec. II D 1.
We next turn to the Fourier transform of S, which will yield an p PSF. Including in its definition a factor (2π) omitted in Eq. (18), we have Here and henceforth, the index +1 is identified with 1.
Since the dependence of K on frequencies enters through the variables ω 1···i , it is convenient to express S through these, too. The conditions ω i = E i+1 i implied by the δ functions are equivalent to ω 1···i = i j=1 E j+1 j = E i1 , and particularly ω 1··· = E +1 1 = E 11 = 0. We thus obtain Our Eqs. (21) and (22) are compact Lehmann representations for the PSFs, with ω i serving as placeholder for E i+1i , and ω 1···i for E i+11 , respectively. We now insert Eqs. (20b) and (22a) into Eq. (18) for G(ω). The corresponding G(ω), extracted as in Eq. (16), has the form of an ( −1)-fold convolution: As already mentioned, for each of G, K, and S, the frequency arguments are understood to sum to zero. Next, we consider the case of an arbitrary permutation p in Eq. (17). The Fourier transforms of functions with permuted arguments, K(ω p ) and S[O p ](ω p ), readily follow from those given above for K(ω) and S[O](ω). For a given p: i → i, the integration measure and exponent in the Fourier integral are invariant under relabeling times and frequencies as t i → t i and ω i → ω i , i.e. d t e iω·t = d t p e iωp·tp . Hence, the above discussion applies unchanged, except for the index relabeling i → i on times, frequencies, and operator superscripts. Depending on context, it may or may not be useful to additionally rename the dummy summation indices as (1 2... ).
For instance, the permuted version of Eq. (21) can be written in either of the following forms: The first choice ensures that the density matrix carries the same index, ρ 1 , irrespective of p. This is helpful for the NRG implementation [34], used to obtain the numerical results in Sec. IV. The second choice yields matrix elements O i ii+1 whose subscript indices are linked to the superscripts. This is often convenient for obtaining analytical results. It shows, e.g., that PSFs whose arguments are cyclically related are proportional to each other. To be explicit, let p = (1... ) and p λ = (λ... 1...λ−1) be cyclically related permutations, e.g. (1234) and (3412). According to Eq. (24b), the corresponding PSFs, S[O p ](ω p ) and S[O p λ ](ω p λ ), differ only in the indices on the density matrix, which appears as ρ 1 or ρ λ , respectively; they otherwise contain the same product i , written in two different, cyclically related orders. Since ρ λ = ρ 1 e −βE λ1 in thermal equilibrium and the δ functions enforce E λ1 = ω 1···λ−1 , we obtain the cyclicity relation This relation is useful for analytical arguments and, in particular, allows one to reduce the number of PSFs in the spectral representation from ! to ( −1)!. However, we here refrain from doing so, since it would increase the complexity of the kernels and modify their role in the spectral representation by introducing Boltzmann factors.
We are now ready to present our final results for the Fourier transform of Eqs. (17). It involves a sum over permuted versions of Eq. (23): Here, ω p = (ω 1 ,..., ω ) is a permuted version of ω = (ω 1 ,..., ω ), with ω 1··· = ω 1··· = 0 understood, and the integral is over the first −1 independent components of ω p . Alternatively, the integration variables can also be chosen independent of p (by using the same set of −1 independent components of ω for all p), in which case the measure will be written as d −1 ω . The permuted ZF convolution kernels [Eq. (20c)] are given by Equations (26) to (28) give the spectral representation for ZF p correlators. Combining them, we obtain For = 2, this reproduces Eq. (4); see App. A. We conclude this section with a remark on connected correlators. These are relevant in many contexts, particularly for the 4p vertices in Sec. IV. The connected correlator G con follows from the full correlator G by subtracting the disconnected part G dis ; the latter is a sum over products of lower-point correlators. Through the spectral representation, we can transfer the notion of a connected and disconnected part from G onto the PSFs, S con = S−S dis . By linearity, G con will follow by combining S con with the same kernels K as for the full correlators.
Let us consider explicitly the case of four fermionic creation or annihilation operators (cf. Sec. III). Then, since expectation values of an odd number of operators vanish. Each factor on the right describes a timedependent 2p PSF, e.g.

C. Matsubara formalism
The derivation of spectral representations in the MF is analogous to that in the ZF, and it utilizes the same real-frequency PSFs from Eq. (28). All arguments regarding time-translational invariance and energy conservation from above still hold. Hence, we also use the same notation in terms of G and G. Nevertheless, there are subtleties arising in the MF. The first concerns (composite) bosonic frequencies that may be zero and lead to anomalous terms as in Eq. (6) [38,39]. The second stems from the fact that imaginary times can be chosen positive and thus always larger than the time set to zero. This leads to only ( −1)! permutations of operators. Yet, the nontrivial boundary conditions of the imaginary-time integral provide terms, thus making up for the total of ! terms. For us, it will be convenient to directly work with ! summands. In fact, it turns out that only the terms arising from the trivial lower boundary of the imaginary-time integral contribute; all others cancel out and need not be computed explicitly. Without this trick, the calculations are more tedious, but still straightforward; see App. B for = 3 and 4.
We start from the MF p function, with operators time ordered on the interval τ i ∈ (0, β) (cf. Fig. 2) and periodic or antiperiodic boundary conditions for G, depending on the bosonic or fermionic nature of the corresponding operators. We wish to compute where the Kronecker δ follows from translational invariance in τ . As for the ZF, we express G(τ ) as a sum over permutations involving products of step functions θ: where Ω p = iω p −ω p , and ω 1··· = 0, set by the δ functions in S, is understood implicitly from Eq. (35b) onward.
where the integrals are arranged in "descending" order, Since Eq. (33) for G(iω) contains a factor βδ ω 1··· ,0 , the integrals in Eq. (37) must yield a result of the form Inserted into Eq. (35b) for G(iω), the K term reproduces the discrete δ function in Eq. (33), yielding while the K term must yield zero when summed over p, The K term in Eq. (38) comes from the outermost integral over τ in Eq. (37), involving e Ω 1··· τ = e iω 1··· τ (recall ω 1··· = 0). However, this integral generates the prefactor βδ ω 1··· ,0 only for terms with no further dependence on τ . These terms, which arise solely from the lower integration boundaries of all subsequent integrals, give K: Conversely, all terms arising from one or more upper integration boundaries, which depend on τ via τ i+1··· , contribute to K . We need not compute them explicitly, since, by Eq. (40), their contributions cancel. We now evaluate the integrals (41) for K. We temporarily exclude the anomalous term arising if some exponents vanish and indicate this by putting a tilde on K (and G): Inserting Eqs. (42) and (28) for S into Eq. (39), we obtain This compact expression is our first main result for MF p functions. For = 2, it yields [upon relabeling summation indices as in Eq. (24b)] For applications, we will be mostly interested in fermionic systems and p functions with ≤ 4. These can still contain bosonic operators, e.g., as bilinears of fermionic operators. However, for arbitrary 2p functions, 3p functions with only one bosonic operator, and fermionic 4p functions, ω 1···i , with i < , produces at most one bosonic frequency. In this case, at most one frequency in the denominators of Eq. (43) can vanish, say ω 1···j , for some j < . One can show (see App. B) that the corresponding anomalous terms are fully included using the following kernel in Eq. (39): Here, δ Ω 1···j ,0 is symbolic notation indicating that the anomalous second term contributes only if Ω 1···j = 0 (i.e. if both ω 1···j = 0 and ω 1···j = 0, the latter due to a degeneracy, E j+1 = E 1 , in the spectrum). In this case, 1/Ω 1···j diverges, but the product i =j duly excludes such factors. The regular first term needs no such exclusion, since, if Ω 1···j → 0, then also Ω i+1··· → 0, and the 1/Ω 1···i divergence is canceled by a −1/Ω i+1··· divergence stemming from a cyclically related permutation. We confirm this in App. B by treating nominally vanishing denominators as infinitesimal and explicitly tracking the cancellation of divergent terms. This procedure yields the following expression for the full kernel: The kernels (45) and (46) yield equivalent results. The former is well suited for analytical work; the latter is more convenient for numerical computations, as it is manifestly free from divergences. The spectral representation given by Eqs. (39), (28), and (45) or (46) is our final result for MF p functions.

D. Keldysh formalism
The Keldysh formalism (KF) [16,17] is based on ordering operators on a doubled time contour involving a forward and a backward branch. In the contour basis, each time argument carries an extra index specifying which branch the corresponding operator resides on. The Keldysh basis employs linear combinations of such contour-ordered correlators. Before discussing these two options in turn, we introduce a fully retarded kernel, a very useful object through which all other KF kernels can be expressed. While doing so, we carefully discuss the imaginary parts of complex frequencies, needed for convergence of real-time integrals. We present a choice of imaginary parts which is consistent both if they are infinitesimal, as assumed in analytical work, or finite, as needed for numerical computations.

Fully retarded kernel
Operators on the forward branch of the Keldysh contour are time ordered while those on the backward branch are anti-time-ordered. In Eq. (17b) of the ZF, we already encountered the time-ordered kernel. In the KF, it will be useful to have a kernel that combines η−1 anti-timeordering and −η time-ordering factors in the form for 1 ≤ η ≤ . As usual, a product over an empty set, with lower limit larger than upper limit, is defined to equal unity. Note how the doubled time contour of the KF (cf. Fig. 4 below) is reflected in Eq. (47): the successive (from right to left) time-ordering and anti-time-ordering factors single out t η , the ηth component of t p , as largest time.
To Fourier transform K [η] (t p ), we first consider the identity permutation, i = i, requiring no subscripts p or overbars. (The general case will follow by suitably reinstating overbars at the end.) In the Fourier integral, we can take the perspective that the largest time t η runs over the entire real axis, while all other t i are constrained by t i < t η . It is thus natural to use the integral over t η for energy conservation and, exploiting time-translational invariance, consider all other time dependencies t i −t η as 1 [4] ω 2 [4] ω 3 [4] ω 4 [4] ω 3 [4] ω 1 [4] ω 4 [4] ω 2 [4] advanced (i.e. contributing only for t i −t η < 0): The t i =η integrals can be regularized, without affecting the t η integral, by replacing the real tuple ω by a complex tuple ω [η] with components ω [η] i , having finite imaginary parts. We thus shift ω i → ω [η] i [see Fig. 3 i . (49) This assigns a negative imaginary part to each frequency ω [η] i =η multiplying t i −t η , as in the 2p case, whereas ω [η] η has a positive imaginary part. The superscript indicates that this choice depends on η. By construction n . For numerics, γ 0 > 0 should be kept small but finite. For analytical arguments, it can be taken infinitesimal.
Since K [η] depends only on time differences, with t η singled out as the largest time, we take t η and t i = t i −t i+1 for i < as our integration variables. Then, having used ω [η] 1··· = 0 in the last step. Now, consider a general permutation p. To compute K [η] (ω p ), the Fourier transform of K [η] (t p ), we need permuted versions of the complex frequencies in Eq. (49). We ) as the tuple (ω [η] ) p obtained by permuting the components of ω [η] , including their imaginary parts, according to p. This moves ω [η] η , the component with positive imaginary part, to slot µ = p −1 (η) [see Fig. 3(b)]. Hence, the components of ω [η] p have a positive imaginary part in the slot i = µ, with the entry i = µ = η, and negative imaginary parts in all other slots: 1··· = 0, and the imaginary part of ω [η] 1···i is negative for 1 ≤ i < µ and positive for µ ≤ i < (yielding To achieve convergent Fourier integrals for K [η] (ω p ), we should thus add negative imaginary parts to all frequencies except the Fourier partner of t η , i.e. ω η , sitting in slot η of ω p . This is achieved by using the complex tuple ω [η] p . Indeed, with superscript η, the positive imaginary part ends up in slot This is our main result for general retarded kernels, applicable for both finite or infinitesimal imaginary parts. Among the denominators, η − 1 ( − η) of them have a negative (positive) imaginary part, for all permutations p. Indeed, for infinitesimal imaginary parts, Eq. (52) reads However, the selection of the components of ω that receive a positive or negative imaginary shift in Eq. (53) depends on p. A specific component ω i may receive a positive shift for one permutation and a negative shift for another. If η = 1, e.g., only denominators of the form ω + 1···i , with positive imaginary parts, occur in Eq. (53). For a permutation with 1 = 1, one of them equals ω + 1 ; with = 1, another equals ω + 1··· −1 = −(ω − 1 ). Hence, a simple sum p K [η] (ω p ) has no well-defined analytical structure w.r.t. ω 1 , or more generally w.r.t. ω. However, we will later obtain suitable combinations of retarded kernels with different η that do, yielding a retarded correlator.
To conclude this section, we note that the above approach to construct K [η] can also be used to find the ZF kernel with finite imaginary parts, as needed for numerics. Indeed, for η = 1, Eq. (47) exactly matches the ZF time-ordered kernel of Eq. (17b), K [1] = K. This kernel has both a largest time t 1 and a smallest time t , and the energy-conservation integral may be performed using either (i) t 1 or (ii) t . Choice (i) corresponds to the above discussion; hence, the η = 1 version of Eq. (52) can be used for K(ω p ) for an expression with finite imaginary parts instead of Eq. (27). For choice (ii), the dependence on all t i −t , i < , is retarded (contributing only for t i −t > 0), as in Eq. (20a). To regularize the corresponding time integrals there, we replace the frequencies ω i in Eqs. (20) by ω [ ] i * , such that those with i < all have positive imaginary parts. The two choices can be summarized as since, for i < , we have γ = −iγ 0 for any permutation. Thus, both choices give positive imaginary shifts accompanying all ω 1···i (i < ), consistent with ω 1···i +i0 + used in Eq. (27). They are both legitimate and yield equivalent results for γ 0 → 0. For = 2, they are identical even at finite γ 0 . For > 2 and finite γ 0 , they lead to qualitatively similar results, with slight quantitatively differences that decrease with γ 0 . For the curves shown in Fig. 6, e.g., both choices yield indistinguishable results.

Contour basis
In the contour basis, time-dependent operators O i (t ci i ) are defined on the Keldysh double time contour. The contour index c i on the time argument t ci i specifies the branch, with c i = − or + for an operator residing on the forward or backward branch, respectively. Correspondingly, a contour-ordered p function, defined as carries contour indices, c = c 1 ··· c , one for each operator.
The contour-ordering operator T c rearranges the operators such that those on the forward branch are all time ordered, those on the backward branch anti-time-ordered; the former are applied first, i.e., they all sit to the right of the latter. It also provides a sign change for each transposition of two fermionic operators incurred while reordering. As for the ZF, we express G c (t) as In this sum over permutations of indices, each summand is a product of a kernel K c , enforcing contour ordering, and the PSF S of Eq. (17c), containing the time-dependent operators. For a given choice of times t and contour indices c, only one permutation yields a nonzero result, namely the one which arranges the operators in contour-ordered fashion. Given c, let λ denote the number of its + entries, i.e. the number of operators on the backward branch. Contour ordering places all λ operators on the backward branch 4. KF time ordering. The -tuples of times t = (t1,..., t ) and contour indices c = (c1,..., c ) specify an -tuple (t c 1 1 ,..., t c ), with t ∓ i on the forward or backward branch, respectively. Consider such an -tuple, with λ contour indices equal to +, the others −. Then, only that permutation p in Eq. (56) yields a nonzero contribution, K cp (tp) = 0, denoted In the corresponding operator product S[Op](tp), all (forward-branch) t − times appear to the right of all (backward-branch) t + times, with larger t − to the left of smaller t − , and smaller t + to the left of larger t + . In other words, the right-to-left order of operators in the product matches the order in which their times appear on the contour.
to the left of all −λ operators on the forward branch. Hence, only those components of K cp (t p ) are nonzero for which c p = c 1 ···c has the form +···+−···−, with λ entries of + followed by −λ entries of −, 0 ≤ λ ≤ . We use the shorthand c p = [λ, −λ] for this structure, e.g. [0, 3] = −−− or [2, 2] = ++−−. Then, the successive time-ordering and anti-time-ordering rules on the two branches (from right to left) are implemented by the kernel cf. Fig. 4. Again, a product over an empty set, with lower limit larger than upper limit, equals unity. The superscript of K, standing for c p = [λ, − λ], reflects the number of + and − entries in c and hence is the same for all nonzero components K cp associated with a given K c .

Keldysh basis
Correlators in the Keldysh basis G k carry Keldysh indices, k = k 1 ···k , with k i ∈ {1, 2}. They are obtained from correlators in the contour basis by applying a linear transformation D to each contour index [17], with c = − or + giving the first or second column of D, respectively. In Eq. (56), the dependence on c resides solely in K c ; thus, the Keldysh rotation yields (with a conventional prefactor) To perform the sum over all c p in Eq. (62b), we recalled that the kernels K cp are nonzero only if c p has the form c p = [λ, −λ], with λ ∈ [0, ]. For these, √ 2D k i c i equals (−1) k i for i ≤ λ and 1 otherwise, yielding the factor (−1) k 1···λ . We used k 1···λ = λ i=1 k i , as usual defining the sum over an empty set as zero, k 1···0 = 0.
Using Eq. (58), we can express the permuted KF kernel through fully retarded ones: This result directly yields the Fourier transform (ω p ) given by Eq. (52). Thus, we know all ingredients of the Fourier transform of G k , written as a convolution of K k and S: Together, Eqs. (64), (63), and (52) give the spectral representation for KF p functions in the Keldysh basis, in a form well suited for numerical computations [34].
To obtain more analytical insight, it is fruitful to further simplify Eq. (63). The fraction in its last line vanishes whenever k λ = 1. Hence, there is a cancellation pattern that depends on the number of 2's, say α, contained in the composite Keldysh index k p = k 1 ···k . To elaborate on this, we first consider the identity permutation, i = i, for which the Keldysh indices of K kp match those of G k . Let η j denote the slot of the jth 2 in k, with η j < η j+1 . Then, k is uniquely specified by the ordered list k = [η 1 ... η α ], e.g. 1111 = [ ], 2121 = [13]. From Eq. (63), it follows immediately that K [ ] = K 1···1 = 0. Next, consider Keldysh indices containing a solitary 2 in slot η, i.e. k η = 2 and k = [η]. As anticipated by this notation, K k is then equal to the retarded kernel K [η] . Indeed, Eq. (63) with k = [η] has only one nonzero summand, having λ = η, and since k i =η = 1 implies k 1···η−1 = η −1. Finally, Keldysh indices with 2 ≤ α ≤ many 2's in slots [η 1 ... η α ] yield a term of the form (65) for each η j : To find the sign, we used k i =ηj = 1, k ηj = 2, and the fact that the η j 's are ordered, implying k 1···ηj −1 = η j +j −2. Now, consider a general permutation p. The permuted Keldysh index k p contains the same number of 2's as k, but in different slots. We can analogously express it through an ordered list as k p = [η 1 ...η α ], whereη j denotes the slot of the jth 2 in k p , withη j <η j+1 . To find theη j 's given the η j 's, note that a 2 from slot η j in k is moved by p to slot µ j = p −1 (η j ) in k p . The sequence [µ 1 ...
In App. C, we show how Eqs. (67) reproduce the wellknown results for 2p KF correlators. Since the set ofη j is obtained by ordering the set of µ j = p −1 (η j ), it follows that eachη j in Eq. (67b) is equal to some η j (with j and j related in a manner depending on p). Hence, the external frequencies ω enter the kernel through sets of complex frequencies ω 1···i , whose imaginary parts are determined by the external Keldysh indices k = [η 1 ... η α ].
The spectral representation (67) constitutes our main result for KF p functions. It is very compact, with the number of terms increasing with α, the number Keldysh indices equal to 2, and offers insight into the analytical structure of Keldysh correlators, as we explain next.
For the correlators G [η] with a solitary Keldysh index 2 in slot η, Eq. (67b) for the kernel K [η] involves only one summand, similar to its analogs in the ZF and MF (without anomalous terms), Eqs. (54a) and (42), respectively. The ZF, MF, and KF results read respectively, in terms of the PSFs S from Eq. (28). There is an important distinction between the two real-frequency correlators of Eqs. (68a) and (68c): for the former, the imaginary part of each component of ω [1] depends on the permutation p; for the latter, the imaginary parts of ω [η] are independent of p, being determined by the external index η for all permutations. Hence, G [η] can be expressed through a single set of complex frequencies ω [η] (as stated before), while G cannot.
Comparing G [η] to the MF correlator (68b), we observe that they exactly agree up to a replacement of frequencies: This is the analytic continuation between MF p functions and (fully) retarded ones in the KF [25,45,46]. It generalizes the well-known 2p relation G 21/12 (ω) = G(iω → ω ± ). By contrast, Keldysh correlators with multiple 2's cannot be obtained from MF ones by direct analytic continuation because their kernels in Eq. (67b) involve two or more sets of frequencies, ω [η1] , ω [η2] , etc., having different imaginary parts. Therefore, they do not have any well-defined regions of analyticity in the space of complex frequencies.
One may also realize that the summation in Eq. (67b) for α > 1 combines denominators of frequencies differing only by their infinitesimal imaginary parts and thus leads to δ functions. These are clearly at odds with any simple analytic structure. Nonetheless, the spectral representations derived above offer a convenient starting point for a systematic analysis of relations between different Keldysh correlators, similar to the famous fluctuation-dissipation theorem for = 2 [cf. Eq. (C4)] [47][48][49]. Such an analysis exceeds the scope of this work and will appear in a separate publication [26].

III. QUANTITIES TO COMPUTE
The next two sections are devoted to illustrating the potential of our approach for computing p functions with several exemplary applications. In the present section, we recall the definition of the 4p vertex, give the Hamiltonians of the relevant models, and summarize various analytical results available for them. In Sec. IV, we present numerical results for the 4p vertex, comparing them against benchmarks where available.

A. Definition of the 4p vertex
For the numerical evaluation of our spectral representations, we focus on local 4p functions of fermionic creation or annihilation operators, d † σ or d σ , where σ ∈ {↑, ↓} is the electronic spin. We also need the 2p correlator, which follows from Eqs. (15), (32) In this section (in contrast to previous ones), we display only the first −1 time or frequency arguments of p functions, evoking time-translational invariance to set the th time argument to zero and ω = −ω 1··· −1 .
In systems with spin symmetry, the 2p function, or propagator, is diagonal in spin and simply reads For the 4p function, in the MF and analogously in the ZF [replacing (−1) 3 by (−i) 3 ]. We here use the Fourier exponent i i ω i τ i with the same sign for all frequencies, whereas one typically attributes alternating signs to annihilation and creation operators. We can directly switch to the standard convention by expressing our final results in the particle-hole representation of frequencies [50] (cf. Fig. 5), involving minus signs for the frequencies ω 2 and ω 4 , related to creation operators. We will use these new variables in the discussion of our results. (An exception is the model for x-ray absorption for which different operators O are used; see Sec. III B for details.) However, for brevity, we retain the ω i notation throughout this section. The 4p correlators can be decomposed into a connected (con) and disconnected (dis) part, G = G con + G dis . The latter corresponds to the independent propagation of two particles and reads in the MF and ZF, respectively, By contrast, the connected part focuses on the mutual interaction. Yet, one can still factor out the propagation of each particle to and from the "scattering event" ("external legs" in a diagrammatic language, see Fig. 5). This yields the effective interaction, i.e. the (full) 4p vertex, where ω 4 = −ω 123 , and analogously in the ZF. Hence, once the disconnected part is removed-which can also be done on the level of the PSFs, see Eq. (30)-the vertex follows from a simple division of numbers. In the KF, each operator is placed on either the forward or backward branch, according to its contour index c i ∈ {±}. Thereby, an p function acquires 2 components. The Keldysh rotation exploits the fact that not all of these components are independent. For = 2, the resulting matrix structure in Keldysh indices k i ∈ {1, 2} is in terms of the retarded (R), advanced (A), and Keldysh (K) component. The matrix structure naturally carries over to Eq. (73). Because of G 11 = 0, the right hand side vanishes if the 4p function has only one Keldysh index equal 2; i.e., the fully retarded components G [η] have no disconnected part. The translation between connected correlator and vertex from Eq. (74) now involves matrix multiplications, i.e., with summation over k i . One thus gets F k σσ from G con;k σσ by multiplying matrix inverses of the propagator (75). With reference to Fig. 5, the right Keldysh index k 1 of G k1k 1 (ω 1 ) corresponds to a creation operator and marks the beginning of the propagator line; the left one, k 1 , corresponds to an annihilation operator and marks the end of the propagator line. Using G 1111 = 0, G 11 = 0 [cf. Eq. (63)] in Eq. (76) directly implies F 2222 = 0. One further finds that a vertex with only one Keldysh index equal 1 in slot η, dubbed F [η] , is directly proportional to G [η] . Hence, we call them (fully) retarded as well. Indeed, in these cases, only retarded or advanced propagators with k i = k i contribute to Eq. (76). When using Eqs. (76) to numerically extract F from G con by dividing out the external legs, the same imaginary frequency shifts must be used for the external-leg 2p correlators on the right as for the 4p correlator G con on the left. In App. D, we explain how this can be achieved.

B. Models
We compute local 4p vertices for three impurity models. The first describes x-ray absorption in metals, the second is the symmetric Anderson impurity model (AIM), and the third is a self-consistent AIM for the one-band Hubbard model (HM) in dynamical mean-field theory (DMFT).
For x-ray absorption, we consider the following Hamiltonian, to be called Mahan impurity model (MIM) [51] (mimicking the nomenclature customary for the AIM): where c = c . The first term describes a conduction band of spinless electrons with flat density of states (1/2D)θ(D−| |) and half-bandwidth D, the second a localized core level with p −D, filled in thermal equilibrium. An x-ray absorption process, described by c † p, transfers an electron from the core level into the conduction band, thereby turning on a local attractive scattering potential −U < 0 with U | p |, described by the third term. We define the absorption threshold ω th as the difference between the ground-state energies of the subspaces with or without a core hole (ω th is of order | p |, but slightly smaller, since the hole-bath interaction is attractive). The x-ray absorption rate at an energy ω relative to the threshold is given by the imaginary part of the particle-hole susceptibility, the 2p correlator χ(ω) = G[p † c, c † p](ω+ω th ). The corresponding 4p correlator is G[p † , c, c † , p](ν, ν , ω), where, in the present context, the definition (72) of these frequencies is replaced by For the arguments equated to ω 1 and ω 4 = −ω 123 , associated with the operators p † and p switching from the no-hole to the one-hole subspace and back, we split off ±ω th , the energy differences E 21 and E 14 associated with the transitions 1|p † |2 and 4|p|1 [cf. Eq. (22b)] between subspace ground states. Furthermore, the bosonic frequency ω is chosen to have opposite sign compared to Eq. (72), ensuring that ω 12 = −ω 34 = ω+ω th matches the argument of the susceptibility G[p † c, c † p](ω + ω th ).

The AIM is described by the Hamiltonian
It contains a band of spinful electrons, a local level with energy d and Coulomb repulsion U , and a hybridization term, fully characterized by the hybridization function ∆(ν) = π|V | 2 δ(ν − ). We take d = −U/2 and choose ∆(ν) = ∆ θ(D − |ν|) box shaped. The local density of states is given by the standard spectral function S std (ν) associated with the 2p correlator

The one-band HM is a lattice model with
where ij enumerate nearest-neighbor lattice sites, t is the hopping amplitude and U the interaction strength.
In DMFT, the HM is mapped onto a self-consistently determined AIM [2]. The associated impurity degrees of freedom d ( †) σ experience the same local interaction U , while their coupling to the rest of the lattice is encoded in the hybridization function ∆(ν). In this paper, we consider the Bethe lattice with infinite coordination number, which yields a semicircular lattice density of states of halfbandwidth D (∝ t), and paramagnetic, spatially uniform phases. Then, the DMFT self-consistency condition simply reads ∆(ν) = (D/2) 2 πS std (ν). Upon self-consistency, local correlators of the HM can be found from the corresponding AIM and thus analyzed in direct analogy.

C. Analytic benchmarks
The MIM was thoroughly investigated by Nozières and collaborators in the ZF [52][53][54]. The particle-hole susceptibility has a power-law divergence for frequencies above the threshold, the celebrated x-ray-edge singularity [51,55], solved analytically [54] as Here, ω is the absorption frequency above the threshold ω th , and δ is the conduction-electron phase shift induced by the core-hole scattering potential. It has the analytic expression [54] δ = arctan(πg), with g = U/(2D) here. It can also be computed with NRG through δ = π∆ h , where ∆ h is the charge drawn in toward the scattering site by the core hole [56,57].
In the first paper of the series [52], the (2p) susceptibility χ was deduced from the (full) 4p vertex. The latter actually contains a variety of power laws, as summarized by Eq. (35) of Ref. 52. Plus, one can extract χ, and thus the same power law (81), from the vertex by sending suitable frequencies to infinity. To this end, we consider the vertex F , related to the 4p correlator G[p † , c, c † , p], whose bare part is F 0 = −U . In the particle-hole representation of frequencies (72), χ then follows as [21,58,59] [60] The limit of fermionic frequencies must be taken such that |ν|, |ν |, and |ν ± ν | become arbitrarily large [21]. In this limit, Im F/U 2 = Im χ(ω) ∝ ω −α . Analytic results are also available for the half filled AIM in the MF in the limits of either weak or infinitely strong interaction. In the former, one often simplifies D → ∞ and T = 0 to find the bare particle-hole susceptibility [61] Its particle-particle counterpart yields −χ 0 . Thus, combining all three two-particle channels, the vertex in secondorder perturbation theory follows as (using↑ = ↓,↓ = ↑) Here, the first term U δ σ,σ on the right is the MF bare vertex F 0;σσ . An expression analogous to Eq. (84) holds for the local vertex of the weakly interacting HM. In that case, Eq. (83) must be computed for the appropriate hybridization function ∆(ν), with a nontrivial frequency dependence resulting from the self-consistency condition.

IV. NUMERICAL VERTEX RESULTS
In the following, we present results obtained by numerically computing local 4p functions with NRG. Generally, NRG allows one to construct a complete basis of approximate eigenstates of the Hamiltonian [67,68] and thus directly evaluate the spectral representations (see Refs. 29 and 30 for computing retarded 2p functions). In the accompanying paper [34], we develop a new NRG scheme to treat 3p and 4p functions. We refer interested readers to Secs. IV-V of that paper for how PSFs are computed with NRG as sums of discrete δ peaks, and to Sec. VI for how these are broadened to smooth curves and how connected correlators and the full vertex are obtained from them. Here, we show the final results of the 4p vertex, for the MIM in the ZF and for the AIM, with both a boxed-shaped and a DMFT self-consistent hybridization, in the MF and KF.

A. ZF for MIM: Power laws
The MIM is a prototypical model for the ZF; see e.g. Refs. 52-54. For this model, NRG proved very successful in computing power-law divergences of 2p functions [27,56,57,67,69,70], but it was never used to investigate similar singularities in 4p functions. As explained in Sec. III, we can extract the famous power law ω −α in terms of the (bosonic) transfer frequency of the particlehole susceptibility χ directly from the 4p vertex F by setting the fermionic frequencies ν, ν to very large values. Figure 6(a) shows various cross sections of Im F (ν, ν , ω), where |ν|, |ν |, and |ν ±ν | are much larger than the halfbandwidth D. All of them collapse onto the same curve.
This curve meets two consistency checks. First, as expected from Eq. (82), it matches Im χ(ω) (red dashed line), computed separately as a 2p correlator. The discrepancy at ω U is due to the broadening of discrete PSFs (described in Ref. 34, Sec. VI.B); it can be removed by reducing broadening, but then wiggly discretization artifacts appear. Second, as expected from Eq. (81), Fig. 6(a) shows a power-law divergence for ω close to the threshold. The exponent α = 0.38 in ω −α (black dashed line, guide to the eye) was obtained from Eq. (81) using the phase shift δ = arctan(πg) = 0.67. We obtained the same value, δ = π∆ h = 0.67, when computing the ∆ h , the charge drawn in to the core hole, with NRG following Ref. 56.
Next, we probe further power laws in the vertex F by setting one fermionic frequency to a large value, ν max , and The real part of F at ω < 0. If only one of ν and ν is large, F follows two distinct scaling behaviors for νmin < −ω < νmax and −ω < νmin, given, respectively, by Eqs. (36) and (40)  the other one to a small value, ν min . As |ω| is reduced from |ω| > ν min to |ω| < ν min , F crosses over between two power laws, given by Eqs. (36) and (40)  These consistency checks, with matching results for highly nontrivial 4p and 2p functions, on the one hand, and agreement between numerical 4p results and analytic predictions, on the other hand, provide confidence that NRG is well suited to compute local 4p functions in the ZF. Moreover, it successfully meets the particularly tough challenge of the regimes ω |ν| |ν | and ω |ν| |ν |, namely resolving a small frequency with exponential accuracy while also keeping track of two larger ones.

B. MF for AIM: Temperature evolution
Most numerical work involving nonperturbative 4p functions is obtained in the MF, where, thanks to the steady progress in quantum Monte Carlo (QMC) techniques, local 4p functions can nowadays be computed with high precision (see [50] for a recent list of references). Using the spectral representation as given by Eqs. (39), (28), (46), they can be computed with NRG, too. For the next parts, we consider the half-filled AIM with box-shaped hybridiza- We start with a moderately low temperature T = 10 −2 D. Figure 7 shows our MF NRG results for the two spin components of F σσ as a function of ν and ν at ω = 0 in the particle-hole representation (72). One can still see the discrete nature of the Matsubara frequencies iν ( ) ∈ iπT (2Z+1). Furthermore, one observes the typical structure in the frequency dependence of the full vertex [20,21,50,63] composed of a background value (independent of ν ( ) ), a distinct signal on the diagonal and antidiagonal (though weaker), and a plus-shaped feature (ν ( ) = ±πT ) (the latter is more pronounced at lower T ). Note that F ↑↑ vanishes identically for ν = ν . This is intuitive from the Pauli principle since, then, all quantum numbers of both fermions involved match. It also follows from the symmetry relation mentioned below Eq. (85), since exchanging either ω 1 ↔ ω 3 or ω 2 ↔ ω 4 at ω = 0 leaves the diagonal (ν = ν ) invariant. At a temperature T = 10 −2 D, we can compare our results to highly accurate QMC data [10,71,72]. We find that the results differ on the level of 1%, which confirms the reliability of our new NRG scheme at moderately low temperatures.
Typical QMC algorithms scale unfavorably with inverse temperature. An important advantage of our MF NRG scheme is hence that it extends to arbitrarily low T . For the given parameters, the Kondo temperature is T K 5×10 −3 D [73]. Accordingly, T = 10 −2 D, as used above, is not low enough to enter the strong-coupling regime, but T = 10 −4 D, used for Fig. 8, is. There, we find that the features already observed in Fig. 7 become sharper and more pronounced. Particularly interesting is the region ν ( ) T K , describing the Fermi liquid with an impurity screened by the Kondo cloud. Indeed, the inset in Fig. 8 shows that the vertex is strongly reduced in this regime, thus giving way to weakly interacting quasiparticles.
To explore this concept further, renormalized perturbation theory (RPT) offers a way of extracting the quasiparticle interaction (local Landau parameter) directly from the NRG low-energy spectrum [75][76][77][78]. This proceeds by comparing the eigenenergies of states with two excited quasiparticles to those with only one; it does not require knowledge of frequency-dependent correlation functions. The resulting valuesŨ σσ =Ũ δ σ,σ should match the lowenergy limit of the effective interaction (i.e. the 4p vertex F σσ with all frequencies sent to zero), multiplied by the quasiparticle weight Z:Ũ σσ = Z 2 F σσ (iω → 0) [76]. Figure 9 shows the vertex evaluated at the lowest Matsubara frequencies, ν ( ) = ±πT , ω = 0, as a function of decreasing temperature. At very large temperatures, T D, correlation effects are suppressed, and the vertex reduces to the bare interaction, F σσ → U δ σ,σ . As we lower temperature much below D, there are strong renormalization effects, and particularly F ↑↓ for ν = ν and F ↑↑ for ν = −ν grow in magnitude. (Recall that F ↑↑ vanishes identically for ν = ν .) Now, for T on the order of the Kondo temperature T K , this trend comes to a halt, and the low-energy components of the vertex start to decrease again. This is the nonperturbative crossover from strongly interacting particles to weakly interacting quasiparticles, as one enters the Fermi-liquid regime. Indeed, for temperatures sufficiently below T K , we find that the low-energy vertex has precisely the same form as the bare vertex: it is nonzero only for different spins, with a value independent of the signs of the frequencies ν ( ) = ±πT . This value precisely matches the RPT estimate, F σσ →Ũ σσ /Z 2 . Vice versa,Ũ σσ can be determined from Z 2 F σσ (iω → 0). Both the decrease of F σσ (iω) for |ω i |, T < T K and the small quasiparticle weight Z 0.36 lead to a small (but finite) quasiparticle interaction ofŨ 0.20U .
To our best knowledge, comparing RPT to the lowenergy limit of a nonperturbative vertex computation has not been realized before. It is a stringent consistency check for the underlying Fermi liquid and confirms the reliability of our NRG scheme down to very low temperatures.

C. KF for AIM: Benchmarks and strong coupling
To our best knowledge, nonperturbative results for KF 4p functions have not been obtained in the literature before. Hence, we begin our analysis with two benchmark FIG. 9. The vertex evaluated at the lowest Matsuabra frequencies ν ( ) = ±πT and ω = 0, denoted by F σσ (±1, ±1, 0), as a function of decreasing temperature. For T D, the vertex reduces to the bare interaction F 0;σσ = U δ σ,σ . Upon lowering T , it exhibits strong renormalization effects and undergoes a crossover from increasing to decreasing magnitude for T > TK and T < TK, respectively. For T TK, it converges to the quasiparticle interactionŨ σσ =Ũ δ σ,σ , divided by twice the quasiparticle weight Z. The parametersŨ 0.20U and Z 0.36 were found from RPT and NRG, independent of the 4p computation, and thus provide a strong consistency check.
cases. The first concerns the limit of infinitely strong interaction, ∆/U = 0, i.e., the atomic limit of the AIM or simply the Hubbard atom (HA). In this case, the many-body basis consists of only four states, and a NRG computation reduces to a simple, exact diagonalization. The PSFs involve only a few δ singularities, which, at half filling, are placed at multiples of U/2. Any real-frequency correlation function directly inherits these singularities. In numerical calculations, one sets a minimal imaginary part γ 0 in the spectral representation. As the denominators contain three factors of the type (ω i +iγ 0 −ω i ), the poles reach a magnitude of γ −3 0 . The MF vertex of the HA is well understood analytically [9], see Eq. (85), and reveals many features shared by general MF vertices; see Fig. 10. We now analyze the KF vertex of the HA. To this end, we compute the spectral representation involving 24 permutations, subtract the disconnected part, and amputate the external legs (cf. Sec. III). As the calculation in this limit is exact, these steps pose no further difficulties. The result are 16 Keldysh components, each having a real and imaginary part. Figure 11 shows the huge variety of features that can be observed as a function of ν and ν at ω = 0 [in the particle-hole representation (72)]. We display five different Keldysh components; all others, related by permutations of their Keldysh indices, show analogous features. A very compact, analytic result for the fully retarded components F [η] can be deduced from the known Matsubara result (85) and the analytic continuation (69). Our numerical results match those to floating-point precision. As explained in Sec. III C, the typical diagonal features of the vertex become δ functions in the atomic limit. They do not appear in F [η] but in other Keldysh components; see Fig. 11. For these other components, too, analytic results have recently been obtained [64], and they perfectly match our numerical ones.
Next, we turn on the coupling to the bath. As we now have a continuous spectrum, discrete PSFs, obtained from a finite number of terms in the sum of Eq. (28), must be broadened. The subtraction of the disconnected part G dis σσ , which requires exact cancellations of large terms, is then numerically difficult. Note that, even if the retarded components G [η] do not have a disconnected part from an analytical perspective, this again relies on exact cancellations that are easily violated numerically. To obtain the vertex via G con σσ most accurately, we employ a twofold strategy [34]. First, we compute G con σσ directly from S con = S − S dis (as discussed at the end of Sec. II B), where S dis is deduced from the full S through appropriate sum rules and subtracted prior to broadening. This eliminates the disconnected part to a large extent but, due to imperfect numerical accuracy, not completely. Second, we use a Keldysh analog of the equation-of-motion method as presented in Ref. 79. It is based on expressing G con σσ through auxiliary correlators, obtained by differentiating G σσ w.r.t., say, the first time argument. Clearly, the choice of this time argument introduces a bias. While this does not appear important in MF applications, we found it to be highly relevant in the KF. As a consequence, the connected part of only those Keldysh components that similarly single out one time argument, i.e. G [η] , could be reliably determined. For other components, we expect a symmetric version of the equation of motion to be beneficial, similar to the one presented in Ref. 80 for the MF. This is, however, left for future work.
Our second KF benchmark involves the limit of weak interaction, U/∆ = 1/2. There, one has an analytic result from second-order perturbation theory (PT), conveniently obtained in the wide-band and zero-temperature limit, The left (right) panels show the real and imaginary parts for σ = σ (σ = σ ). The top row shows NRG data, with U/D = 1/20 and T = 10 −3 D; the bottom row shows results of second-order perturbation theory (PT) in the wideband and zero-temperature limit. We find very good qualitative agreement (the overall magnitude depends on the broadening prescription), even though the limit of weak interaction is highly nontrivial for the diagonalization-based NRG. which can be found by analytic continuation of Eq. (84). One should keep in mind that weak interaction, and even the noninteracting limit, is highly nontrivial for a diagonalization-based algorithm like NRG. In Fig. 12, we compare our NRG results of the weakly interacting AIM to PT, considering the real and imaginary parts of F [1] σσ , as a function of ν and ν at ω = 0 for both spin components. While the precise values of the NRG vertex depend on the broadening prescription, we overall find very good qualitative agreement between NRG and PT. There are dominant diagonal structures that arise whenever one of the (bosonic) frequency arguments in (the analytic continuation of) Eq. (84) vanishes. According to Eq. (83), their real part describes a peak of width ∼ ∆, as reproduced in Fig. 12. Upon subtraction of the bare vertex, the ↑↓ component has no background contribution, since χ 0 (ω 12 = ω) in Eq. (84) comes with a factor δ σ,σ . The ↑↑ component has a background value of χ 0 (0); its diagonal vanishes upon cancellation of χ 0 (ω 12 = ω) and χ 0 (ω 14 = ν−ν ). Whereas these features are easily explained from a perturbative perspective on the vertex, obtaining them in a numerical approach that starts from the spectral representation of the correlator is a stringent test for the whole machinery.
After passing both of these benchmarks, we can confidently present our results for the retarded vertex in the strongly interacting regime. We use identical parameters as above: U/∆ = 5, U/D = 1/5, and T = 10 −4 D below the Kondo temperature T K 5×10 −3 D. In Fig. 13(a), we show Re F [1] ↑↓ −F [1] 0;↑↓ as a function of ν and ν , for three choices of ω. At large ω, the vertex mostly involves only one diagonal. This structure can again be understood from a diagrammatic perspective [21] related to Eq. (84): contributions from both the particle-hole channel corresponding to χ 0 (ω 12 = ω) and the particle-particle channel, corresponding to χ 0 (ω 13 = ω +ν +ν ), are suppressed at large ω, while those of the other particle-hole channel persist independent of it, similar to χ 0 (ω 14 = ν −ν ). As we decrease ω, more features build up: Re F [1] ↑↓ develops an increasing background value, independent of ν and ν , except for distinct values forming a plus shape. Figure 13(b) shows the real and imaginary parts of F [1] σσ − F [1] 0;σσ , for both spin components, at ω = 0. We again enlarge the low-energy window ν ( ) T K (dashed square). Similar to the MF results (cf. Fig. 9), the lowenergy limit of F [η] should reproduce the RPT prediction. Figure 13(c) shows the difference F [1] σσ −Ũ σσ /(2Z 2 ) (the factor of 2 comes from the Keldysh rotation). Indeed, this difference is an order of magnitude smaller than F [1] σσ −F [1] 0;σσ , and, for ν ( ) → 0, it goes to zero. While the real parts of both spin components show a rather extended regime of small values, the imaginary parts, particularly regarding ↑↓, exhibit only a thin line of zero values. 1D |ν|, |ν | 0.5D, similar to where the standard (2p) spectral function S std (ω) has dips separating quasiparticle peak and Hubbard bands (Fig. 14). The insulating solution has very sharp diagonal lines (ν = ν ).

D. MF and KF for HM: DMFT solution
For the strongly interacting AIM, the MF vertex in Fig. 8 and the real part of the retarded KF vertex in Fig. 13 look somewhat similar. [Recall that KF 4p functions inherit a factor 1/2 from the Keldysh rotation; see e.g. Eq. (69).] This is drastically different for the results we present in the following. There, we analyze the (half filled) one-band HM within DMFT, which amounts to solving a self-consistently determined AIM [2]. We take U/D = 2.6 and T /D = 10 −4 in the coexistence region of a metallic and insulating solution. Figure 14 shows the standard (2p) spectral function S std of the strongly interacting AIM [ Fig. 14(a)], considered previously, and the metallic and insulating solution of the self-consistent AIM describing the HM [ Fig. 14(b)]. The self-consistent metallic solution has much more pronounced features, where the spectral weight between the quasiparticle peak and the Hubbard bands almost goes to zero. The insulating solution has a gap at ω = 0 and broad Hubbard bands around the positions of the atomic peaks, ±U/2. Figure 15 shows the local 4p vertex for the DMFT solution of the HM, displaying the MF vertex on the left, and a retarded Keldysh component, k = 1222, on the right. In the top row, depicting the metallic solution, the MF vertex is reminiscent of the AIM results in Fig. 8, albeit with significantly larger values and a more extended plus-shaped structure. By contrast, the KF vertex reveals entirely new features not found in the AIM results of Fig. 13. Next to the typical structure consisting of a background value, diagonal line, and plusshaped structure, it exhibits sign changes at intermediate frequencies 0.1 |ν ( ) |/D 0.5. They are thus at similar frequencies as the aforementioned dips in S std (Fig. 14) and likely related to these. In the bottom row, depicting the insulating solution, the MF vertex is very similar to the HA solution (Fig. 10), with almost divergent values. The KF vertex also looks similar to that, but different from the KF HA results in Fig. 11. The reason is that, for the HA, the diagonal structure expected in F [1] has a width that vanishes as ∆ → 0 and enters as a δ function in other Keldysh components. However, for the insulating solution of the HM, the hybridization, though gapped, remains finite. It thus allows for a diagonal signal in F [1] , even if it is only a very sharp line.
The overall magnitudes of the KF results in Fig. 11 are rather different from the MF results. On the one hand, this may result from the fact that the real-frequency results do have a richer structure in this case. On the other hand, the precise values of the KF results still depend on the broadening prescription, particularly for the insulating solution. Further improving the broadening of 4p real-frequency data is planned for follow-up work. This will enable a thorough analysis of the DMFT realfrequency vertex and promises valuable insight into strongcorrelation effects on the two-particle level.

A. Summary
The many-body problem is typically addressed by either deriving -point ( p) correlation and vertex functions from an action or by working with operators and states in relation to a Hamiltonian. The connection between both approaches through spectral or Lehmann representations is well known for various 2p functions and for imaginary-time functions with ≤ 4. Here, we completed this picture in a generalized framework, providing spectral representations for arbitrary p functions in three commonly used many-body frameworks: the real-frequency zero-temperature formalism (ZF), imaginary-frequency Matsubara formalism (MF), and real-frequency Keldysh formalism (KF).
Through the spectral representations, we elucidated how p correlators G in the ZF, MF, and KF are related to one another. We expressed G as a convolution of partial spectral functions (PSFs) S and kernel functions K. The PSFs are formalism independent and contain the dynamical information of a given system. By contrast, the kernels are system independent and encode the time-ordering prescriptions of the formalism at hand. We first derived the spectral representation in the ZF where it is most compact. We proceeded with the MF, using analogous arguments and the same PSFs, and discussed anomalous terms that arise when vanishing (bosonic) Matsubara frequencies and degeneracies in the spectrum occur together. For the KF in both the contour and Keldysh bases, we identified a (fully) retarded kernel K [η] through which all other KF kernels can be expressed. Among the KF correlators, those with a solitary Keldysh index equal to 2, G [η] , have the simplest spectral representation. It precisely matches the one of the MF (without anomalous parts) up to a replacement of imaginary frequencies by real frequencies with infinitesimal imaginary parts, iω → ω [η] , making the analytic continuation between MF and retarded KF p functions manifest.
We used a novel NRG method, described in the accompanying paper [34], to evaluate the spectral representations of 4p functions for selected quantum impurity models. Starting with a simple model for x-ray absorption treated in the ZF, we analyzed multiple power laws in the 4p vertex. Proceeding with the Anderson impurity model (AIM) in the MF, we benchmarked our NRG results against Monte Carlo data at intermediate temperatures.
The NRG technique is beneficial at very low temperatures, where we found the vertex to exhibit strongly reduced values for Matsubara frequencies below the Kondo temperature T K . Indeed, by decreasing the temperature T from above half-bandwidth D to below T K , we studied the crossover from strongly interacting particles to weakly interacting quasiparticles on the level of the 4p vertex F (effective interaction). At T T K , we deduced from F the quasiparticle interaction, which has a form identical to the bare electronic interaction, with values predictable from renormalized perturbation theory (RPT) and the NRG low-energy spectrum.
For the AIM in the KF, we first tested our NRG results in the solvable limits of weak and infinitely strong interaction. Finding qualitative agreement with perturbation theory for weak interaction was a demanding test for the diagonalization-based NRG. For infinitely strong interaction, i.e. the Hubbard atom (HA), we found rich features in the different Keldysh components of the vertex. In the intermediate regime of the strongly interacting AIM, we observed the formation of renormalization effects in the retarded vertex F [1] with decreasing the transfer frequency ω and showed that F [1] reproduces RPT for temperature and frequencies below T K . Overall, the retarded vertex showed features that might be expected from the MF result and perturbation theory. This was drastically different in our final results where we analyzed the one-band Hubbard model within DMFT in the regime of coexisting metallic and insulating solutions. For the strongly correlated metal, the KF vertex showed distinct features at intermediate frequencies, absent in the MF counterpart. We speculate that these go hand in hand with the dips separating the quasiparticle peak and the Hubbard bands in the standard (2p) spectral function. Regarding the insulating solution, the MF vertex is similar to the MF HA result. The KF vertex also looks similar to that, but the finite hybridization, albeit gapped, leads to notable deviations from the KF HA.

B. Outlook
In this work, we have focused on (i) the formal properties of spectral representations and (ii) benchmark tests of our NRG scheme [34] to compute the 4p vertex of quantum impurity models. This sets the stage for a variety of intriguing applications.
Let us start with more formal aspects. In the Introduction, we mentioned the divergences of two-particle irreducible vertices in the MF. Mathematically, these divergences arise through the matrix inversion of a generalized susceptibility χ, a function of discrete Matsubara frequencies, whose eigenvalues become negative and thus cross zero at some point [10]. Expressed through real frequencies, the former sums over matrix elements become integrals over complex functions. At finite temperature, there is an additional matrix structure of Keldysh indices. The tools presented here allow one to compute χ in the ZF or KF, and to investigate if divergences of the 2PI vertex persist for real frequencies.
Moreover, real frequencies provide the natural language for Fermi-liquid theory. Indeed, the original Fermi-liquid works used the zero-temperature approach, defining the Landau parameters in terms of the ZF (full) 4p vertex [82]. In the mean time, the T = 0 Landau parameters were expressed through Matsubara vertices, too [11,83,84]. Such MF vertices can also be computed at T > 0. However, a stringent extension of Fermi-liquid theory to finite temperature should ideally use real-frequency Keldysh objects. Their nonperturbative evaluation becomes accessible through this work.
On the practical side, an important application of local 4p functions is given by diagrammatic extensions of DMFT. We briefly recall that DMFT describes local correlations, assuming a purely local self-energy. Diagrammatic extensions of DMFT employ diagrammatic relations to further incorporate correlations of arbitrary wavelength [50]. For instance, in the ladder dynamical vertex approximation [42,85,86], dual fermion formalism [87][88][89], or a functional renormalization group (fRG) flow [90] starting from DMFT [91], a momentum-dependent self-energy and vertex are constructed from the local (full) 4p vertex. Many successful results along these lines have been obtained in the MF [50]. Yet, the commonly used Monte Carlo-based DMFT impurity solvers are not able to reach very low temperatures. Our MF NRG results, which extend to arbitrarily low temperature, can remedy this limitation. Furthermore, to properly interpret the intriguing phases of strongly correlated electron systems, real-frequency diagrammatic extensions of DMFT would be invaluable. The momentum dependence could be generated by real-frequency diagrammatic techniques, such as the rather well-developed Keldysh fRG [65,92,93]. The necessary building blocks are real-frequency local 4p functions, such as the KF NRG results shown here.
Another interesting topic requiring the computation of real-frequency 4p functions is the theory of resonant inelastic x-ray scattering (RIXS) of strongly correlated materials [94]. We present some proof-of-principle RIXS results in the accompanying paper [34].

ACKNOWLEDGMENTS
We thank Anxiang Ge, Gabriel Kotliar, and Andreas Weichselbaum for insightful discussions. We also thank A. Ge for making his analytical KF results for the Hubbard atom available to us as benchmark data. Furthermore, we thank Alessandro Toschi and Patrick Chalupa for fruitful discussions and for providing the QMC benchmark data. Discussions with Daniel Springer and Clemens Watzenböck in early stages of this work are appreciated. For our numerical computations, we exploited symmetries using the QSpace library developed by A. Weichselbaum [31,33] For brevity, we wrote G(ω) = G(ω), K(ω p ) = K(ω 1 ), and S(ω p ) = S(ω 1 ), hiding the second frequency argument, with ω 1 = −ω 2 = ω and ω 2 = −ω 1 implicit. The ingredients needed above are K(ω i ) = 1/ω + i and, using Eq. (24b), Inserting Eqs. (A2) into Eq. (A1), we obtain Eq. (4): In the main text, we derived the MF formulas somewhat indirectly, arguing that all contributions from the upper integration boundaries in Eq. (37) cancel. For completeness, we here give a direct derivation of the MF Lehmann representations for = 2, 3, 4. We also discuss anomalous terms arising when denominators vanish. For this purpose, we focus on correlators for which at most one frequency ω 1···i , with i < , is bosonic, and derive the expressions (45) and (46) for the full MF kernel K(Ω p ) given in the main text. The reasons for this focus are stated before Eq. (45); other cases can be treated analogously.
In the following, we consider the cases = 2, 3, and 4 in turn. For each, we start by computing the regular contributions, signified by a tilde on k and G, for which all denominators arising from the τ i integrals are assumed to be nonzero. We subsequently discuss anomalous cases, signified by a hat on k, for which this assumption does not hold. For these, we recompute the corresponding integrals more carefully. Alternatively and more elegantly, the anomalous terms can also be found from the regular ones by a limiting procedure, treating nominally vanishing denominators as infinitesimal rather than zero.
Since the divergences from the denominators of the regular k cancel after summing over permutations [cf. Eq. (B9)], we may also represent the full kernel as corresponding to Eq. (45). Here, δ Ω 1 ,0 is symbolic notation indicating that the anomalous second term is nonzero only if Ω 1 = 0. No restriction is placed on the Ω 1 values for the first term, with the understanding that vanishing denominators should be treated using infinitesimals. (B13) Regular part: Performing the integral (B13) and assuming all resulting denominators to be nonzero, we obtain k(Ω (123) ) = 1 For the last line, we combined the β-independent terms using (1/Ω 1 − 1/Ω 12 )/Ω 2 = 1/(Ω 1 Ω 12 ). We also exploited Ω 123 = 0 to rewrite the denominators of the β-dependent terms in a way that reveals the denominators of all three terms to be cyclically related by (312) → (231) → (123).
Since ω 12 = 0 implies S (312) = S (123) , we can express the product k(Ω (123) )S (123) as − Here, we split the anomalous contribution equally between the first and third terms, using Ω 3 = 0 to rewrite the denominators such that their Ω's match the nonzero Ω's in the first and third terms of Eq. (B14b) for k. Jointly, Eqs. (B19) and (B14b) imply that G(iω) of Eq. (B2) has the form (39), with a kernel defined as if Ω 1 Ω 12 = 0, for any of the six permutations p = (123). This is consistent with the general Eq. (46) from the main text.
Since the anomalous terms not proportional to β stem from the numerators of the regular part k, they need not be displayed separately-they are generated automatically when treating vanishing denominators as infinitesimal and summing over permutations. The full kernel can thus be expressed in the following form, equivalent to Eq. (B20), but written in the notation of Eq. (45):  The sum over all six permutations q = (1234) yields G(iω) = p ζ p dω 1 dω 2 dω 3 S (1234) where the sum now runs over all 24 permutations p = (1234). This is again consistent with Eqs. (39) and (42).
Anomalous part: If all four operators are fermionic, ω i and Ω i are always nonzero, but Ω ij can vanish. Indeed, in Eq. (B23), denominators vanish (i) in the second and fourth terms if Ω 12 = −Ω 34 = 0, and (ii) in the first and third terms if Ω 23 = −Ω 41 = 0. We discuss these cases by taking the vanishing frequencies to be infinitesimal; recomputing the integral (B22) yields the same results.
Finally, we remark that there are attempts in the literature to incorporate anomalous contributions known from the MF, as the one in χ 0 (iω), into retarded functions like χ R 0 (ω) (see Ref. 39, and references therein). This is not necessary when working within the full-fledged KF, as the information of anomalous MF contributions is encoded in other Keldysh components, such as χ K 0 (ω). Indeed, the second-order retarded self-energy, yields the correct result Σ R (ν) = U 2 /(4ν + ) upon using χ A 0 (ω) = [χ R 0 (ω)] * = 0 and χ K 0 (ω) = 2πiδ(ω). Any artificial modification of χ R 0 would give an incorrect result.