Computing local multipoint correlators using the numerical renormalization group

Local three- and four-point correlators yield important insight into strongly correlated systems and have many applications. However, the nonperturbative, accurate computation of multipoint correlators is challenging, particularly in the real-frequency domain for systems at low temperatures. In the accompanying paper, we introduce generalized spectral representations for multipoint correlators. Here, we develop a numerical renormalization group (NRG) approach, capable of efficiently evaluating these spectral representations, to compute local three- and four-point correlators of quantum impurity models. The key objects in our scheme are partial spectral functions, encoding the system's dynamical information. Their computation via NRG allows us to simultaneously resolve various multiparticle excitations down to the lowest energies. By subsequently convolving the partial spectral functions with appropriate kernels, we obtain multipoint correlators in the imaginary-frequency Matsubara, the real-frequency zero-temperature, and the real-frequency Keldysh formalisms. We present exemplary results for the connected four-point correlators of the Anderson impurity model, and for resonant inelastic x-ray scattering (RIXS) spectra of related impurity models. Our method can treat temperatures and frequencies -- imaginary or real -- of all magnitudes, from large to arbitrarily small ones.


I. INTRODUCTION
Correlation functions beyond the standard two-point (2p) functions, such as 3p and 4p functions, play a central role in many-body theory. For instance, the 4p vertex describes the effective interaction between two particles in a many-body environment and signals pairing instabilities; 3p functions give the full detail of a particle reacting to an external perturbation and connect to 2p correlators through various Ward identities. Through their physical significance and mutual relations, 2p, 3p, and 4p functions form the basis of the microscopic Fermi-liquid theory.
In the accompanying paper [1], we describe a novel approach for analyzing and computing -point ( p) correlators. Its central idea is to separate system-dependent spectral information from formalism-dependent analytic properties. This is achieved via the convolution of partial spectral functions (PSFs) with various types of convolution kernels. We consider three different frameworks, the imaginary-frequency Matsubara formalism (MF), the realfrequency zero-temperature formalism (ZF), and the realfrequency finite-temperature Keldysh formalism (KF). The convolution kernels are responsible for the analytic structure of the theory; they depend on the choice of formalism, but are system independent. The PSFs encode the spectral properties of the system; they are independent of the choice of formalism and can be viewed as the central porters of information about the dynamics of the system. They have generalized Lehmann representations, expressed through the eigenenergies and corresponding many-body eigenstates of the system. Having compact support, bounded by the largest energy scale in the system, they are particularly convenient for numerical storage.
The present paper offers a detailed discussion of how NRG can be utilized to compute 3p and 4p correlators for quantum impurity models. To this end, we generalize the full-density-matrix (fdm) NRG approach [21,22], devised for 2p correlators, to the 3p and 4p cases. To tackle the novel dependence of p PSFs on −1 frequencies, we introduce an additional, iterative scheme to finely resolve regimes involving frequencies of different magnitudes, |ω i | |ω j | (1 ≤ i, j < ). The bulk of this paper is devoted to a detailed description of this approach.
We have performed numerous benchmark checks of our new method, focusing on the computation of the local connected 4p correlator and corresponding 4p vertex of various impurity models. We obtained excellent agreement with (i) analytical predictions for the power-law behavior of the ZF vertex for x-ray absorption, (ii) QMC results for the MF vertex of the Anderson impurity model (AIM) at intermediate temperatures, (iii) exact results for the KF vertex of the AIM with infinitely strong interactions (the Hubbard atom), (iv) perturbative results for the KF vertex of the AIM with weak interactions, and (v) perturbative results for the ZF and KF connected correlator for the AIM with weak interactions. The numerical results for checks (i)-(iv) were discussed in Ref. 1. Here, we present the underlying PSFs and the results for check (v). The success of all these tests establishes the reliability of our method, and its ability to treat temperatures and frequencies-imaginary or real-of all magnitudes, including very small ones. For example, for the AIM, our approach can reach temperatures much lower than the Kondo temperature and accurately resolve correspondingly low frequencies.
In this paper, we also include an application of great physical interest, the computation of resonant inelastic x-ray scattering (RIXS) spectra. RIXS is a powerful experimental technique for probing various excitations in solids over a wide energy range [53,54]. However, the effects of many-body correlations on RIXS spectra are poorly understood, even for simple models. Numerical calculations of RIXS spectra are typically based either on ED of small systems [55][56][57][58][59][60][61][62][63][64][65][66] or on a Bethe-Salpeter approach building on ab initio calculations using density functional theory [67][68][69]. These methods have limited ability for capturing strong-correlation phenomena characterized by low energies and long length scales, such as the Kondo effect, or the emergence of a small quasiparticle coherence energy scale in many correlated metals.
Our method is ideally suited for overcoming this limitation. We demonstrate this with proof-of-principle calculations of the RIXS spectra of two minimal models: the Mahan impurity model (MIM) involving a free conduction band interacting with a core hole, used to describe x-ray absorption in metals in a seminal work by Mahan [70], and an augmented AIM (AAIM), involving the AIM and a core hole. We elucidate how the RIXS spectra of these models are affected by Anderson orthogonality and the Kondo effect, respectively. Note that the scope of the AIM goes beyond impurity models, as DMFT and its diagrammatic extensions describe lattice systems by the AIM with a self-consistently determined bath.
As an overview, we here summarize the workflow of our NRG method for the case of local 4p correlators. First, one constructs a Wilson chain by discretizing the impurity model and obtains a complete basis of (approximate) energy eigenstates of the entire chain; Sec. III outlines this step. Then, one computes 4p PSFs in a recursive way, described in Sec. V, such that contributions to 4p PSFs involving frequencies of widely different magnitudes can be obtained by invoking the routines for 3p and 2p PSFs. To enable such a recursive approach, we rephrase the established scheme for computing 2p PSFs, introducing additional notational and diagrammatic conventions, as elaborated in Sec. IV.
By convolving the 4p PSFs with the kernels given in Sec. II, one obtains the full 4p correlators in the MF, ZF, and KF. For the real-frequency ZF or KF correlators, it is necessary to broaden the 4p PSFs, which are discrete due to the discretization in the beginning. Further, to describe genuine two-particle properties, such as 4p vertices, one needs to extract the connected part of 4p correlators. This can be numerically challenging. In Sec. VI, we describe strategies for improving the accuracy of these steps.
Our result for connected 4p correlators are presented in Sec. VII, and for RIXS spectra in Sec. VIII. Finally, Sec. IX offers a summary and an outlook to applications opened up by our NRG approach to 3p and 4p correlators.

II. SPECTRAL REPRESENTATIONS
We begin with a summary of the spectral representations of p correlators derived in Ref. 1, defining p PSFs and their convolution kernels. To describe key ideas and introduce notation, we will first focus on the ZF, stating analogous MF and KF results thereafter.

A. Zero-temperature formalism
In the ZF time domain, an p correlator is defined by a time-ordered product, where the kernel K(t p ) = is nonzero only if the permuted times t p = (t 1 , ..., t ) satisfy t i > t i+1 . Here, p (12. with ω = (ω 1 , ... , ω ). Using the shorthand notation ω 1···i = i j=1 ω j , the δ function on the right expresses energy conservation, ω 1··· = 0, following from timetranslational invariance for G(t). Because of the multiplicative structure of Eq. (1), G(ω) can be expressed as an ( −1)-fold convolution of the Fourier transforms of K and the operator product, resulting in an expression of the form Here, G, K, and S each have only − 1 independent arguments, with ω and ω p = (ω 1 ,..., ω ) understood to obey energy conservation, ω 1··· = 0 and ω 1··· = 0. The ZF convolution kernel K can be chosen as with the shorthand ξ i···j = j m=i ξ m . This corresponds to Eq. (54a) of Ref. 1. For numerical calculations, the imaginary parts have small but finite values γ i > 0. The PSFs S have a Lehmann representation: Here, each underlined summation index i enumerates a complete set of many-body eigenstates |i of the Hamiltonian H, with eigenenergies E i , transition energies E ji = E j − E i , and matrix elements O ij = i|O|j , ρ 1 = e −βE1 /Z. (We use calligraphic symbols for operators, roman ones for matrix elements.) Equations (3)- (5) give the spectral representation for ZF p correlators.

B. Matsubara formalism
In the MF, operators time-evolve as O i (τ ) = e Hτ O i e −Hτ , and p correlators are defined as G(τ ) = (−1) −1 T i=1 O i (τ i ) . Operators are time-ordered on the imaginary-time interval τ i ∈ (0, β). The Fourier transform of G(τ ) takes the form with ω denoting a set of discrete Matsubara frequencies, and the Kronecker symbol enforcing energy conservation. A permutation expansion for G(τ ) analogous to Eq. (1) leads to the spectral representation, with real frequencies ω p , and ω 1··· = 0 and ω 1··· = 0 understood. The PSFs S are again given by Eq. (5). The MF kernel, expressed through Ω p = iω p − ω p = (Ω 1 , ..., Ω ) and Ω 1···i = i j=1 Ω j , reads The first line gives the regular part of the kernel, applicable if all denominators are nonzero. If this is not the case, anomalous contributions arise. The second line gives their form under the assumption that at most one denominator vanishes, say, Ω 1···j , for some j < . This includes the case of arbitrary 2p correlators, 3p correlators with only one bosonic operator, as well as 4p correlators of fermionic operators, such that ω 1···i , with i < , produces at most one bosonic frequency. Equations (5), (7), and (8) give the spectral representation for MF p correlators.

C. Keldysh formalism
KF correlators in the contour basis, G c , are defined as where t ci i are real times and T c denotes contour ordering. They carry a tuple of contour indices c = c 1 ··· c , with c i = − or + if operator O i resides on the forward or backward branch of the Keldysh contour, respectively. In this work, we will treat KF correlators in the Keldysh basis, G k , carrying a tuple k = k 1 ··· k of Keldysh indices k i ∈ {1, 2}. The correlator G k is obtained from G c by applying a linear transformation D to each contour index, The Fourier transform of G k is defined as in Eq. (2). The resulting G k (ω) has the spectral representation, involving real frequencies ω and ω p , with ω 1··· = 0 and ω 1··· = 0 understood. The PSFs S are yet again given by Eq. (5). The KF kernel can be expressed as

D. Structure of partial spectral functions
A very attractive feature of the above spectral representations is that all three formalisms, ZF, MF, and KF, contain the same PSFs, given by Eq. (5). These PSFs are thus the central porters of dynamical information. The main goal of this paper is to describe a numerical algorithm for computing them. We conclude this section with some general remarks on their structure.
First, for a given permutation p, the kernels K and the PSFs S all depend on the −1 integration variables (ω 1 , ... , ω −1 ) only in the combinations ε i = ω 1···i , for 1 ≤ i < . For numerical purposes, it is convenient to use these combinations as independent integration variables. We collect them in an ( −1)-tuple, ε = (ε 1 , ... , ε −1 ) and express the PSFs through this tuple, defining We will henceforth display the subscript on S p (ε) only in formulas involving a permutation sum, such as Eq. (9); elsewhere the subscript p will be suppressed, it being understood that S(ε) refers to a specified permutation. When computing S(ε), the spectral resolution can be controlled individually for each ε i . This is crucial for NRG computations, where the spectral resolution attainable for each ε i is not uniform but proportional to |ε i |. For a given permutation p, the index contraction pattern of S(ε), as given by Eq. (5) with ε i = ω 1···i there, can be depicted diagrammatically as follows (for = 4): Here, the matrix elements O i ii+1 and ρ 1 are represented by black squares and a dot, respectively. An arrow i on a solid line connecting two such symbols denotes a sum over i. It points away from the second (ket) index and toward the first (bra) index of the corresponding operator matrix elements. Each δ function is represented by a dashed line, whose arrow points from the incoming index 1 to the outgoing index i + 1 in the corresponding condition ε i = E i+1 − E 1 . It will be convenient to express S in a form reflecting the nested structure of the dashed lines: Here, O ε is defined as an operator with matrix elements with ( ij ) ε a shorthand notation for multiplication by a δ function. Similarly, the matrix elements of ( O 1 ) ε1 O 2 ε2 involve a combination of δ functions with nested arguments: This nesting complicates the numerical computation of S using NRG, as will be elaborated in the next sections.

III. WILSON CHAINS
In this section, we summarize the key ingredients of Wilson's NRG approach for solving quantum impurity models. We begin with the construction of Wilson shells, sets of approximate energy eigenstates resolving successively lower-energy parts of the spectrum with ever finer resolution. We then review how complete sets of eigenstates for the entire system can be constructed from the Wilson shell eigenstates. Finally, we give a brief preview of how these complete sets of states can be used in Lehmann representations to compute a 2p PSF. This sets the stage for a more detailed discussion of the computation of PSFs for = 2, 3, and 4 in subsequent sections.

A. Wilson shells
A quantum impurity model describes a discrete degree of freedom coupled to a continuous bath of excitations. In the NRG approach devised by Wilson [19], the bath is discretized on a logarithmic grid characterized by a discretization parameter Λ (> 1). The discretized bath levels have energies typically scaling as ∼ ±Λ −k D, with integer k ≥ 0 and D the half-bandwidth (taken as unit of energy). Then, the model is mapped onto a semiinfinite "Wilson chain," where the impurity is represented by the leftmost site, having index −1, and the bath by the sites 0, 1,... . Site n has on-site energy n and couples to site n−1 via hopping amplitude t n . The exponential decay of t n ∼ Λ −n/2 and the even stronger decay of n ensures energy-scale separation. In practice, one studies finite chains, choosing the largest site index N sufficiently large to resolve all relevant energy scales of the problem, including the temperature (requiring Λ −N/2 T ). The Hilbert space of such a "length-N " chain is spanned by the local basis states {|σ N } = {|σ N ⊗ · · · ⊗ |σ −1 }, with local dimension d per site for sites n ≥ 0. (Following Ref. 22, we build the direct product basis from right to left, as if filling the chain in the order c † N ... c † −1 |0 , where for brevity we displayed only site indices.) Exploiting energy-scale separation, the chain can be diagonalized iteratively, starting with the impurity and increasing the subchain length n, one site at a time. This amounts to resolving ever smaller energy scales, resulting in a topdown (high-to low-energy) refinement strategy. However, the subchain Hilbert space dimension, ∼ d n , increases exponentially. Hence, a truncation scheme is needed, in which some states are kept (K) and others discarded (D).
Let n 0 be the last site for which the Hamiltonian of a length-n 0 subchain H n0 can be diagonalized exactly without truncation. The set of eigenstates of H n0 is known as Wilson shell n 0 . Its lowest-lying levels have characteristic level spacing ∼ Λ −n0/2 , as that is the smallest energy scale in H n0 . Now, one partitions this shell into two sectors, containing low-lying kept and high-lying discarded states, respectively, {|s n0 X }, X ∈ {K, D}, with eigenenergies E n0 Xs . Then, one proceeds iteratively as follows.
For any n > n 0 , suppose that H n−1 has been diagonalized, yielding the eigenstates {|s n−1 X }, X ∈ {K, D} of shell n−1, with eigenenergies E n−1 Xs . Now, add site n and diagonalize H n on the truncated space spanned by the direct product of the new site and the kept sector of shell n−1, {|σ n ⊗|s n−1 K }. (The neglect of the high-lying discarded sector during this step is justified by energy-scale separation [19].) The resulting eigenstates form shell n. It has low-lying level spacings ∼ Λ −n/2 ; hence, it provides a refined description of the K sector of shell n−1, which had larger spacings ∼ Λ −(n−1)/2 . If n < N , partition the eigenstates of shell n again into low-lying kept and highlying discarded sectors, {|s n X }, with eigenenergies E n Xs , concluding the nth iteration step. At the last iteration n = N , one declares all eigenstates of H N as discarded, since there is no next iteration in need of kept states.
The eigenstates obtained in step n can be written as linear combinations of the form [22] |s n X = |σ n ⊗ |s n−1 where sums over the repeated indices s, σ n are implied. The coefficients are arranged as matrices, M σn KX , with elements labeled by ss . The matrices M σn and M σn † appearing in the ket (15a) and bra (15b) can be graphically represented by three-leg vertices, with incoming or outgoing arrows associated with the kets on the right or left of Eq. (15a), respectively, and all arrows inverted for the bra version, Eq. (15b). Equations (15) also apply for the first few iterations up to n 0 , just without truncation. (For the impurity site, M σ−1 K is a vector, not a matrix.) The eigenstates of H n thus obtained are matrix product states (MPSs) [71,72] of the form They form an orthonormal set, In our MPS diagrams, sums over s or σ indices, indicated by bonds between vertices, will usually not be displayed but understood to be implicit. Vertices on rightor left-pointing lines for kets or bras represent M 's or M † 's, respectively, lines connecting vertices denote index contractions, and the empty circle indicates an identity matrix. We will henceforth suppress arrows on vertical lines representing local state spaces.

B. Anders-Schiller basis for Wilson chain
Following Anders and Schiller [73], we next construct states defined not on subchains but the full length-N chain. For this purpose, we change our perspective: whereas, above, H n was defined on a length-n chain, we henceforth define it as H n = H N | tn>n= n>n =0 , living on the full length-N chain just as H N , but with couplings t n = 0 and on-site energies n = 0 turned off for all n > n. A topdown refinement of the chain, starting from H n0 , then successively replaces H n−1 by H n , "switching on" the parameters t n , n . By diagonalizing each H n and combining its eigenstates with the discarded states of previous iterations, one obtains a sequence of sets of basis states, each spanning the Hilbert space of the full chain, but with ever finer energy resolution for low-lying energies.
For any subchain ending at site n (≥ n 0 ), we denote the set of states spanning its environment (the rest of the chain) by {|e n = |σ N , . . . , σ n+1 }. Then, the states resolve the spectrum of H n , with low-lying level spacings ∼ Λ −n/2 , but carry an additional degeneracy d N −n due to the environmental degrees of freedom. For shell n 0 , the set of all states {|se n0 X } forms a complete basis for the full length-N chain, albeit with a coarse energy resolution and very high degeneracy. For any later shell n, the union of all its states with the discarded states from previous iterations, ∪ n0≤n<n {|se n D } ∪ {|se n X }, likewise forms a complete basis for the length-N chain. [Completeness is proven below; see Eq. (28).] With increasing n, the low-energy spectral resolution becomes exponentially finer and the environmental degeneracy exponentially smaller. For the last shell at site N , all states are defined as D and there are no environmental states and no environmental degeneracy. The set ∪ n0≤n≤N {|se n D }, comprising the discarded states of all shells, is again complete on the full length-N chain and known as Anders-Schiller (AS) basis.
Evoking energy-scale separation, we now adopt the NRG approximation. It states that, when acting on states from shell n, the Hamiltonian of the full chain, H = H N , can be approximated by that of the first n sites, H n , yielding the eigenvalue E n Xs : The AS basis, collecting all discarded states, thus forms a complete set of approximate eigenstates of the full H, with spectrum {E n Ds } (and degeneracy d N −n for each E n Ds ).

C. Spectral functions from NRG: Brief preview
Here, we briefly preview the computation of 2p PSF by NRG, to set the stage for the formal developments of subsequent sections. Since the AS basis is complete and equipped with eigenenergy information (albeit approximate), it can in principle be used to compute spectral functions via Lehmann representations [21,22], with the identification |i = |se n D . For = 2, e.g., Eq. (12) yields with approximate eigenenergies E n Ds , E n Ds . Since the energy resolution of the AS basis becomes finer with increasing n, the spectral resolution attainable becomes exponentially fine with decreasing energy ε. This fact makes NRG a singularly powerful tool for studying the low-energy dynamics of quantum impurity models.
Yet, the Lehmann representation (20) is not the final expression used for NRG calculations. The reason is that the evaluation of shell-off-diagonal (n = n) contributions to the double sum is computationally demanding, without improving the accuracy of the results. A shell-off-diagonal contribution, say with transition energy ε = E n Ds −E n Ds for n < n, involves the difference between two energy values with different resolutions, ∼ Λ −n/2 and ∼ Λ −n/2 . The frequency resolution of such off-diagonal contributions is dominated by the coarser Λ −n/2 of the earlier shell. The better resolution of the later shell thus is futile, yielding no added benefit-the later shell is overrefined.
We thus arrive at a central principle for the NRG computation of spectral functions: avoid shell-off-diagonal contributions and find shell-diagonal representations. Reserving a detailed discussion for later sections, we here just state the main idea: the K sector of a given shell n may be viewed as a coarse-grained description of all later shells n > n (after all, the latter are obtained by refining the former). Thus, the overrefined off-diagonal contributions to Lehmann representations can be coarse grained bottomup (low-to high-energy) by replacing n>n,se |se n D n D se| by se |se n K n K se|, i.e., the projector onto all shells later than n by the projector onto the KK sector of shell n. As shown in Refs. 21 and 22 and recapitulated in Sec. V A, this leads to a coarse-grained version of Eq. (20) with just a single sum over shell-diagonal contributions: This is the final expression used in NRG calculations.
Here, each transition energy E n Xs −E n Xs involves two energies from the same shell n, computed with the same accuracy ∼ Λ −n/2 . For each n, the sum over sectors includes only DD, DK, and KD matrix elements; KK contributions are excluded, since these are represented, in refined fashion, by later shells. Finally, we note that the frequency dependence of the spectral function is resolved with a resolution comparable to |ε|, furnished by those shells in the expansion for which Λ −n/2 matches |ε|.
The above ideas form the basis of the fdm-NRG approach of Ref. 22 for computing 2p PSFs. Our purpose here is to generalize it to 3p PSFs and 4p PSFs. To this end, we will need various formal properties of Wilson shell projectors, expressed in a compact and economical notation. The next few sections thus develop a formalism for expanding operator products in terms of operator projections to specific shells.
In the material that follows, the expansions of nonlocal operators (Sec. IV C) and their products (Sec. IV F), the operator slicing scheme needed for such expansions (Sec. V B), and the subsequent methods for computing multipoint PSFs and correlators (Secs. V C-VI) are novel methodological developments from this paper. To explain the new notions efficiently, below we will reformulate some well-established ideas and methods, such as the expansion of local operators and the binning [21,22,71], in terms of our formalism.

IV. OPERATOR EXPANSIONS
This section describes how to expand various types of operators and operator products along a Wilson chain. We begin by defining projectors onto specific Wilson shells and discussing their properties. We then consider the expansions of local operators, acting nontrivially only on sites −1 to n 0 . Thereafter, we turn to their time-or frequency-dependent versions, which are nonlocal, in that they depend on the Hamiltonian of the full chain. Next, we discuss the representation of the density operator on the Wilson chain. Finally, we consider expansions of products of two or more operators.

A. Shell projectors
We begin by noting that AS states defined on shells n and n have overlaps of the following forms [cf. Eq. (Here, δ ee refers to the environmental modes of sites > max(n, n).) In words, n = n overlaps vanish unless the "earlier" sector is K, and n = n overlaps yield identity [by Eq. (17)]. These properties will be used repeatedly.
We next introduce shell projectors as convenient tools for working with the AS basis. For any n ≥ n 0 , let denote the projector onto the X sector (with X ∈ {K, D}) of shell n. Since the K sector of shell n spans both the K and D sectors of shell n+1, their projectors satisfy a refinement identity: This can be used to iteratively refine the description of the K sector of shell n in terms of the K and D sectors of subsequent shells up to any shell n, with n < n ≤ N : If n = N , the last term is absent, since we define all states of shell N as D, so that it has no K sector: This identity reflects the fact that the K sector of shell n, spanned by the basis {|se n K }, is also spanned by ∪ n>n {|se n D }, providing a basis for the union of the D sectors of all later shells. The latter basis has a finer energy resolution than the former. Thus, Eq. (26) can be used in two ways: top down, where the right-hand side refines the left, or bottom up, where the left-hand side coarse grains the right. Depending on context, we will adopt either one or the other point of view below. Now, consider the identity operator on the full N -site chain, expressed via the K and D projectors for site n 0 : This resolves the identity using eigenstates of H n0 , but since n 0 typically is small, the energy resolution is poor.
To improve the resolution, we may refine the K sector of shell n 0 via Eq. (24), obtaining a sequence of equivalent resolutions of the identity, for any n from n 0 to N : This proves completeness of the AS basis, and of its less finely resolved analogues, mentioned after Eq. (18). Importantly, any of these bases may thus be used for evaluating the PSFs, as indicated in Sec. III C.
We henceforth view sites n 0 to N as default summation range for n, writing n , n≤ n , and n>n for N n=n0 , n n=n0 , and N n=n+1 , respectively. We will often encounter products of projectors. Such products can be simplified using the following identity, P n X P n X = δ n<n δ XK P n X + δ nn δ XX P n X + δ n>n P n X δ KX , (29) obtained via Eq. (22). The δ symbols indicate that the first, second, and third terms contribute only for n < n, n = n, and n > n, respectively: The projector product identity (29) implies the following rules: a projector product • with n = n vanishes unless the "earlier" sector is K, and then it equals the projector of the later shell; • with n = n is sector diagonal; • with KD or DK vanishes unless K is "earlier" than D; • with DD is site diagonal. These rules will be used extensively in later sections.

B. Expanding local operators
A local operator, acting nontrivially only on sites −1 to n 0 , can be represented exactly in the basis of shell n 0 , The operator O can be expanded in various ways along the Wilson chain. Such expansions involve operator projections to specified combinations of shells and sectors: One obvious expansion strategy sandwiches O between two identities, resolved through the AS basis using Eq. (28): This expansion involves shell-off-diagonal terms. However, for reasons explained in Sec. III C, we prefer an expansion that involves only shell-diagonal terms.
One way to obtain a shell-diagonal expansion is bottomup coarse graining. We replace n>n P n D in Eq. (32) by P n K , using Eq. (26): The expansions (32) and (33) These relations are graphically represented as To be consistent with Eq. (18), we depict the first (second) index of [O n XX ] ss by an incoming (outgoing) leg. We now discuss a more economical, top-down approach, leading directly to the shell-diagonal expansion (33) without a shell-off-diagonal detour. Starting from shell n 0 , we iteratively refine only the KK sector of each successive shell. For a given shell n (≥ n 0 ), we may use Eq. (25) to refine its KK sector in terms of shell n+1, Iterating up to any shell n larger than n, we obtain By construction, all terms generated in this manner are shell diagonal, and each sum on XX always excludes the case KK, except for the site n where the expansion terminates. If it terminates at n = N , where there is only a DD sector, there is no KK contribution at all: where n may take any value from n 0 to N . For the latter choice, shown on the right, Eq. (37) matches Eq. (33). Thus, the top-down and bottom-up approaches for expanding operators are equivalent. The shell-diagonal expansion of local operators is exact, as follows directly from the completeness of the AS basis. For nonlocal operators, the expansion becomes approximate, but it is a reasonable approximation given logarithmic energy resolution ∼ Λ −n/2 . In the next section, we will demonstrate the expansion of nonlocal operators in the top-down approach.

C. Expanding nonlocal operators
Next, we consider time-dependent operators and their Fourier transforms, needed for the computation of spectral functions. Even if O is local, O(t) = e iHt Oe −iHt is not, since its definition involves the full Hamiltonian acting on all sites, H N . The same is true for its Fourier transform O ε = dt 2π e iεt O(t), with matrix elements given by Eq. (13), needed in Eq. (12) for PSFs. Hence, such operators must a priori be expressed through Wilsonchain expansions involving the entire chain. For reasons explained in Sec. III C above, this should be done in a manner leading to a shell-diagonal expansion.
We derive the shell-diagonal expansion of frequencydependent nonlocal operator O ε in the top-down way.
[The expansion of time-dependent operators O(t) can be obtained by Fourier transforming the expansion of O ε .] Hence, we start from a coarse-grained representation of where the Hamiltonian is represented by H n0 . We then refine it by successively turning on the couplings t n and energies n along the chain, i.e., replacing H n−1 by H n . This amounts to repeating the iterative refinement of the KK sector of shell n via Eqs. (35) and (36a), but now with added frequency labels ε: This expansion is approximate as the frequency-dependent matrix elements employ the NRG approximation: Diagrammatically, we represent the δ function in Eq. (39a) by a dashed line bracketing the symbol for the matrix element. The arrow on the dashed line points from the incoming to the outgoing energies in the corresponding condition ε = E n Xs − E n Xs : For a given ε, nonzero matrix elements will come mainly from those shells whose low-lying level spacing ∼ Λ −n/2 is of order ε. However, energy differences between highlying levels of earlier shells will also contribute a bit, since the level spacing within a shell rapidly decreases with increasing energy.
In Eq. (38), the KK projector on the left of each → can be viewed in two ways: as a coarse-grained placeholder, to be replaced by the more refined expression on its right once the need for further refinement arises, or, if there is no such need, as the final expression terminating the expansion. We will henceforth adopt this refinement strategy throughout, but for simplicity use the "=" notation for both local and nonlocal operators. Since the frequency dependence of O ε enters simply via a factor multiplying its matrix elements, we will often drop the superscript, writing just O. Where needed, the frequency dependence can easily be restored.

D. Density matrix
We now discuss the thermal density matrix , a nonlocal operator, needed for thermal averages. In the AS basis, Further employing the NRG approximation, is fully diagonal, where the weight contributed by each discarded state of shell n is given by a Boltzmann factor: The sector projections of for shell n, defined as n XX = P n X P n X and found via Eq. (29), are given by Thus, n XX is sector diagonal. Note that n KK encompasses a sum over all later shells; this sum provides a refinement of shell n in the spirit of Eq. (36b). Hence, for any n from n 0 to N , we can express the density matrix as [cf. Eq.
Next, we consider the reduced density matrix, n XX , for a subchain of length n, obtained from n XX in bottom-up fashion by tracing out the environmental states of all later sites, starting from N and working back to n+1: Equation (42) implies n KD = n DK = 0. We find the matrix elements of n DD and n KK using Eqs. (42) and (44) and the following diagrams, which depict them as circled dots: The vertically aligned dashed lines, understood to be connected, represent the trace over sites > n. For n DD , this trace yields a degeneracy factor d N −n : On the right, Z D n = s e −βE n Ds is the partition function of the D sector of shell n (without its environmental states), and w n = d N −n Z D n /Z the relative weight contributed by that sector to the total partition function. Finally, for n KK , the sum over n > n in Eq. (42) yields Starting at n = N −1, it can be computed recursively via a backward sweep (see above diagram, last line, right-hand part): The reduced density matrix n KK residing in the K sector is used when the density matrix contracts with the shelldiagonal expansion of operators [cf. Eqs. (33) and (38)]. If the projection of an operator onto a K sector does not need further refinement (e.g., since the operator is local), the density matrix to contract with the projection does not need to resolve later shells. Then, the reduced density matrix instead of the original can be used for the contraction through the K sector. We will see an example of this in the next section.

E. Thermal averages
The expectation value O = Tr[ O] of a local operator, nontrivial only on sites ≤ n 0 , can be computed as follows: In the first step, we expressed O through their sector projectors for iteration n 0 , recalling that n0 is sector diagonal and the trace cyclic. Since O is local, it acts as the identity operator on the Hilbert space of all sites > n 0 .
In the second step, we exploited this to trace out these sites, reducing n0 XX to n0 XX . Performing the remaining trace over the Hilbert space of sites ≤ n 0 yields a simple trace over matrix elements (denoted by tr) of shell n 0 , with the eigenstates of shell n 0 made explicit on the right. We will mostly use the more compact trace-of-matrixproduct notation of the middle expression, suppressing s indices. The diagram depicts the expressions on the right of Eq. (46) and Eq. (47), respectively: For future reference, we note that the above strategy also works for thermal averages of shell-n projections: For the first step, we used Eq. (43), = n≤n n DD + n KK , and noted that O n is orthogonal to the D sectors from all earlier shells [cf. Eq. (29)]. The second step mimics Eq. (47), with δ XX from the cyclicity of the trace. Importantly, the above strategy also works for averages of shell-diagonal products, A n B n , A n B n C n , or A n B n C n D n . Such products, just like O n , are orthogonal to all earlier D sectors and trivial on all later shells. Hence, Eq. (48) applies again, with O n XX on the right replaced by a corresponding product of blocks of matrix elements, e.g., XX ρ n XX A n XX B n XX . We will exploit this fact when computing spectral functions below.

F. Expanding operator products
The computation of p PSFs inevitably involves products of frequency-dependent operators, e.g., These require Wilson-chain expansions. Moreover, the latter should be shell diagonal, to avoid overrefinement and to facilitate computing expectation values via Eq. (48). We now explain how such expansions can be obtained.
One approach starts by individually expressing each operator in the product in terms of its full-chain expansion (37) from site n 0 all the way to site N . However, the resulting expansion is not shell diagonal. It can be brought into shell-diagonal form by bottom-up coarse graining the D sectors of shell-off-diagonal terms, but this requires tedious rearrangements, discussed in Appendix A.
A simpler, top-down approach is to iteratively refine entire products rather than individual operators, in a manner that retains shell diagonality throughout. For brevity, we use notation appropriate for local operators.
[For nonlocal ones, frequency labels should be added and the "=" signs in Eqs. (49) and (50) below read as " →" refinements, as discussed at the end of Sec. IV C.] Starting from n 0 , we refine, for each shell n, only the all-K sector of the product, representing it by a sum over all sectors of shell n+1. Using iteration steps analogous to Eq. (35), one finds the following generalizations of expansion (37) for a single operator: We depict shell-diagonal operator products as follows: As an important example, consider the composite op- Operator products as in Eqs. (50) always appear in thermal averages, such that the trace equates the first and last sector index. For = 2, e.g., Eq. (50a) and (48) yield where KK contributions are excluded. However, for products of ≥ 3 operators, the expansions obtained via Eqs. (50) do contain terms involving KK sectors. In the top-down approach, this is an inevitable consequence of insisting on shell diagonality. (In the bottom-up approach, such KK contributions are not present a priori. Yet, they do arise once one coarse grains in order to avoid shell-offdiagonal terms; see Appendix A.) If the corresponding operators are frequency dependent, (O n KK ) ε , their matrix elements [cf. Eq. (39a)] involve δ functions of the form δ(ε − E n KK ), containing KK transition energies. These require further refinement in case ε is much smaller than the characteristic scale of shell n. Strategies for achieving this will be discussed in detail in Secs. V C and V D.
In Eqs. (50), each product is shell diagonal; hence, the energy differences in the δ functions of a given product are all computed with the same accuracy. To be specific, the matrix elements of (A n ) ε1 B n or ((A n ) ε1 B n ) ε2 are given by [cf. Eq. (14)] and analogously for higher products. We depict the resulting contraction patterns as follows: The underlined ε frequencies and slashes through some dashed lines indicate "slicing", a numerical strategy for dealing with products of δ functions discussed in Sec. V B.

V. PARTIAL SPECTRAL FUNCTIONS
We now have all ingredients needed for computing PSFs. We start with the case = 2, recapitulating the fdm-NRG approach of Ref. 22. We then discuss some numerical techniques used to deal with the δ functions arising in Lehmann representations. Finally, we consider the cases = 3 and = 4. We will denote 2p PSFs by S(ε), 3p PSFs by S(ε) with ε = (ε 1 , ε 2 ), and 4p PSFs by S(ε) with ε = (ε 1 , ε 2 , ε 3 ), using the same symbols S throughout. The number of independent frequency arguments, −1, will always be clear from the context.

A. Partial spectral functions: = 2
We begin with the case = 2. By Eq. (12), the 2p PSF of A and B is given by We denote it by S(ε) for short. The above two forms of writing the trace are equivalent and can be used interchangeably. They implement the δ functions occurring in Eq. (5) in different ways while exploiting the fact that is diagonal in the energy eigenbasis. Analogous equivalent forms also exist for the shell-diagonal traces encountered below. In the present section, we focus on the first form; in subsequent sections, we will refer to both. We start from Eq. (55a). To finely resolve the ε dependence of the product ( A) ε B, we use Eq. (50a) to express it as a sum over all shells ≥ n 0 , obtaining the expansion where each S n is a trace over a shell-diagonal product, The following diagram schematically depicts the iterative refinement of KK sectors (red squares) leading to Eq. (56).
The traces in Eq. (57) can be computed using Eq. (48), replacing an operator trace over the entire Wilson chain by a shell-n trace over products of shell-n matrices, and by the reduced density matrix ρ n of shell n: The depiction of Eq. (58) mimics the diagram for Eq. (47): Combining Eqs. (56)- (58) and making the matrix trace over shell n from the latter explicit, we obtain our final formula for 2p PSFs: This formula, restating Eq. (21) in more compact notation, reproduces an expression first given in Ref. 22.
In numerical practice, the Wilson chain is diagonalized in a forward sweep, yielding all coefficients [M n KX ] ss and energies E n sX from site n = −1 to N . Subsequently, the matrix elements of ρ n XX are computed recursively during a backward sweep from site N to n 0 . To evaluate PSFs, the matrix elements A n XX , B n XX and (ρA) n XX are evaluated in a second forward sweep up to site N . The δ function is dealt with by "binning", a technique described next.

B. Binning and slicing
We briefly interrupt our formal development with an interlude on numerical matters. As illustrated by Eq. (59), the Lehmann representations of PSFs involve sums over very many discrete δ functions, originating from the frequency-dependent matrix elements in Eq. (39a). To obtain smooth functions, the discrete δ peaks have to be broadened at the end of the computation through convolution with a suitable broadening kernel, as further explained in Sec. VI B and Appendix C. Therefore, it is not necessary to keep track of the precise position, E n XX,ss , of each peak. (Storing all pairs of spectral weights and peak positions becomes intractable when large numbers of states are kept during the iterative diagonalization.) Instead, we may adopt a binning strategy: we partition the ε axis into narrow intervals such that each interval I ε is centered around ε. Then, for all spectral peaks with positions E n XX,ss lying within I ε , we adjust their positions to ε, i.e., replace δ(ε − E n XX,ss ) by δ(ε − ε). This discretization of frequency variables enables us to store the PSF as a histogram array; each of its elements is the sum of the spectral weights whose transition energies fall into I ε and are then adjusted to ε. This binning strategy is graphically explained in Fig. 1(a).
Because of the logarithmic bath discretization underlying the Wilson chain, the ε grid should have a logarithmic structure, too, fine enough to resolve separate shells. Enumerating grid points ε[m] by an index m, we choose ε[±m] = ±10 (|m|−1)/n dec ε min for m = 0 and ε[m = 0] = 0, where n dec is the number of grid points per decade and ε min = ε[m = 1] is the smallest energy resolved on the grid. We limit the grid size by taking |ε[m]| ≤ ε max , with ε max larger than all other energy scales in the system. The grid point at ε = 0 collects spectral weights sitting at truly zero frequency (up to numerical precision) rather than small but finite frequencies. It is needed for the anomalous MF kernels in Eq. (8) and for constructing the PSFs of disconnected correlators (see Sec. VI A). When the grid is plotted on a logarithmic axis, the spacing between grid points on this axis, ∆ ε = (ln 10)/n dec , should be much narrower than that between the discretized bath energy levels, ∆ LD = ln Λ. In this work we use Λ = 4 and n dec = 32, yielding ∆ ε = 0.072 and ∆ LD = 1.39.
Whereas 2p PSFs have only one frequency argument, the spectral contributions to p PSFs for > 2 have multiple frequency positions, ε 1 , ε 2 ,.... Generally, the frequency values associated with the same spectral contribution can largely differ in size, e.g., |ε 1 | |ε 2 |. In such a case, the value of ε 1 is determined at an earlier shell, that of ε 2 at a later shell. (See Secs. V C and V D for details.) Thus, information on large frequencies associated with early shells has to be recorded while information on smaller frequencies from later shells is still being evaluated. For this, we use a slicing strategy.
For a given shell n and sector XX, the frequencydependent matrix elements ([O n XX ] ss ) ε constitute a threedimensional object labeled by two discrete indices ss and a continuous index ε. Upon discretizing the latter, i.e., replacing ε by ε, this object should be replaced by a rank-three tensor with discrete indices s, s, ε, defined as Each ε value defines a two-dimensional slice through this tensor, containing nonzero entries only if E n XX,ss lies within the bin I ε . Thus, slicing involves no physically motivated approximations-it is merely a convenient way of dealing with a three-dimensional data structure.
The sliced tensor, depicted schematically in Fig. 1(b), has to be stored explicitly. By keeping the ε leg open while computing tensor contractions involving ss, this task can be performed in parallel for all values of ε.
To indicate diagrammatically that a frequency ε is treated by slicing, we will use a perpendicular slash, "slicing" through the associated dashed line, as illustrated in the diagram after Eq. (54). Dashed lines without slashes signify binning. Whether a δ function requires binning or slicing depends on the context; we will explain this in detail in the subsequent two sections.
C. Partial spectral functions: = 3 We next turn to 3p PSFs, S[A, B, C](ε), with two frequency arguments, ε = (ε 1 , ε 2 ). Their computation starts along the same lines as that for = 2, with a top-down shell-diagonal expansion of the threefold operator product using Eq. (50b). However, extra refinement efforts are required if the two frequencies differ significantly, |ε 1 | |ε 2 | or |ε 2 | |ε 1 |. We will express the terms describing these cases through 2p PSFs of a single operator and a composite operator, the latter incorporating the dependence on the larger frequency via slicing. A high resolution for the smaller frequency can then be achieved by a further top-down expansion of the 2p PSFs. This strategy, not needed for = 2, is the main new ingredient for the NRG computation of higher-order PSFs.
We again begin from Eq. (12), or S(ε) for short. The above two ways of writing the trace are equivalent and will be used interchangeably. Expanding Eq. (61b) via Eqs. (50b), we obtain where each S n denotes a shell-diagonal trace: The ε 1 dependence enters via a KK sector and should be further refined. By contrast, the ε 2 dependence, entering via a DK sector, need not be refined further and can be sliced. We indicate this by replacing ε 2 by ε 2 . We thus view B n KD (C n DK ) −ε 2 as the KK sector of a composite operator with a sliced, not direct, frequency dependence. Evoking Eq. (55a), we can then read Eq. (64a) as the ε 1 -dependent 2p PSF of a single and a composite operator: The " " sign signifies that the right-hand side involves the numerical approximation of slicing, which limits the resolution attainable for ε 2 to Λ −n/2 . Note that slicing is needed only for KD and DK sectors, but not for the DD sector, since the latter is excluded during further refinement. This is numerically convenient, since slicing the DD sector would require much more memory than for the KD and DK sectors.
The KDK term can be treated similarly. Starting from Eq. (61a), we express it as ε2 as a composite operator. Its dependence on ε 2 enters through KK elements; hence, we refine it by identifying a 2p PSF via Eq. (55a): We note that there also are other, equivalent ways of associating ε 1 and ε 2 with operators. For example, if S n KDK is expressed via Eq. (61b), then the representation (65b) would be obtained via Eq. (55b). When implementing this scheme numerically, it is convenient to prioritize associating the original frequencies ε i with operators rather than their sign-flipped versions −ε i , unless the sign flip is inevitable, e.g., as in Eq. (64b). This prioritization yields a better organization when the routine for 3p PSFs is invoked recursively, as needed in the next section.
Since the 2p PSFs (64b) and (65b) are built from two KK operators, they can be refined top down to improve the ε 1 or ε 2 resolution, respectively. Iterative use of Eq. (49a) yields expansions similar to Eq. (56): These expansions are restricted to shells n > n, since they are seeded by operators defined on the KK sector of shell n (their XX = KK matrix elements vanish for all shells ≤ n). The expansions (66a) and (66b) yield high resolution of the smaller frequency in the regimes |ε 1 | < |ε 2 | or |ε 2 | < |ε 1 |, respectively. The diagram below depicts the iterative refinement (66b) of the KDK sector for fixed ε 1 , represented by the slanted blue plane (viewed from somewhat below, again with the E 1 -E 2 plane at the front): The middle and right-hand diagrams involve projections of KK composite operators from shell n onto later shells n, e.g., n XX P n X (A n KD ) ε 1 B n DK P n It mimics the diagram after Eq. (34) for finding O n from O n0 , but is seeded by (A n KD ) ε 1 and B n DK , the former sliced. The triple line on the right edge of the effective matrix element signifies a contraction over the ladder from sites n to n. The frequency dependence enters shell diagonally for both ε 1 on site n and ε 2 on site n, with resolution ∼ Λ −n/2 or Λ −n/2 , respectively. For n much larger than n, such matrix elements contribute only for frequencies satisfying |ε 1 | |ε 2 |. In physical terms, they give the amplitudes for the composite operator AB to cause low-energy transitions with energy cost ε 2 via virtual intermediate states with energy cost ε 1 . The ε 2 dependence may be either binned (as here) or sliced (as in Sec. V D).
The expansion (62), with summands S n XX X [A, B, C](ε) given by Eq. (63) for DX X and KDD, Eq. (66a) for KKD, and Eq. (66b) for KDK, constitute our final formulas for S(ε). They express it through n sums of both adequately refined three-operator traces and 2p PSFs enabling further refinement. The numerical cost for computing all adequately refined terms scales as ∼ N . The cost for computing the refinements of all 2p PSFs scales as ∼ N 2 , requiring a n>n sum for every n, but this can be done in a nested fashion, reusing code written for 2p PSFs. In other words, though Eqs. (66) and their grapical depictions are instructive for understanding the structure of the refinement procedure, they need not be coded explicitly; instead, S n KKD and S n KDK can be computed using the 2p PSF subroutine implementing Eq. (59), with some sliced matrix elements as input. Finally, we consider 4p PSFs, S(ε) with ε = (ε 1 , ε 2 , ε 3 ). Evoking Eq. (12), we will use the five equivalent forms There is some freedom in choosing which form to use in a given situation. Our choices below are convenient for expressing some contributions to 4p PSFs through 2p PSFs and 3p PSFs involving composite operators containing sliced matrix elements, facilitating their numerical computation.
Expanding the third form via Eqs. (50c) yields with shell-diagonal traces S n defined as Out of the 15 contributions from the sum over sectors, nine, namely DX X X and KDDD, do not involve products of KK sectors of the three frequency-dependent operators ( A) ε1 , B ε2 , D −ε3 . These nine are adequately resolved and can be computed as shell-n traces using Eq. (48) while binning all three frequencies. By contrast, the other six contributions need KK top-down refinements via PSFs of composite operators. The three terms involving two K's and two D's can be expressed through 2p PSFs of a single operator and a composite triple, or of two composite doubles. For the KKDD term, e.g., we express Eq. (69) through the fourth form in Eq. (67) and then slice the dependence on ε 2 and ε 3 , entering via KD or DK sectors: Here, we view B n KD C n DD (D n DK ) −ε 3 −ε 2 as the KK sector of a composite triple, sliced with respect to both ε 2 and ε 3 , and the second line as a 2p PSF of two shell-n operators [cf. Eq. (55a)]. The KDKD and KDDK terms can be treated analogously, using the second and first forms of Eq. (67): Each of the three 2p PSFs defined in Eqs. (70) depends on a single, nonsliced frequency argument, with the other two frequencies entering via slicing of composite operators. The dependence on the nonsliced frequency can be further refined via a subsequent n>n S n XX expansion, as discussed in earlier sections. The representations used for Eqs. (69) and (70) can be depicted as follows, The remaining three contributions to Eq. (68), involving three K's and one D, can be expressed through 3p PSFs of two single operators and one composite double. Using the fourth, fifth, or first forms of Eq. (67), we obtain identifying 3p PSF by evoking Eq. (61b) for the first two cases and Eq. (61a) for the last one. For each 3p PSF defined in Eqs. (71), the dependence on the two nonsliced frequencies can be refined top down via an expansion n>n S n XX X as in Eq. (62), but restricted to n > n. The resulting terms can be depicted as follows: Each such S n XX X is treated as discussed in Sec. V C: the DX X and KDD contributions via shell-n traces, the KKD and KDK contributions via identification of 2p PSFs and their top-down expansions over shells n > n, thereby achieving independent resolutions for the smaller two frequencies.
The n > n expansion of S n KKKD , e.g., contributes mainly for |ε 3 | > |ε 1 |, |ε 2 |. Each Sñ KKD and Sñ KDK term in this expansion can be refined via a further top-down 2p PSF expansion n>n S n XX , which mainly contributes for |ε 3 | > |ε 2 | > |ε 1 | or |ε 3 | > |ε 1 | > |ε 2 |, respectively: Similarly, the refinement of S n KKDK yields proper resolution for |ε 2 | > |ε 1 |, |ε 3 |, and of S n KDKK for |ε 1 | > |ε 2 |, |ε 3 |. The matrix elements of composite triples arising in such n > n > n expansions are built from nested combinations of matrix elements of composite doubles, such as , and for all KKKD, KKDK and KDKK as ∼ N 3 (a 3p PSF for every n). Though costly, these computations can be performed in a systematic, nested fashion, reusing the subroutines for 2p PSFs when computing 3p PSFs, and the 2p and 3p PSF subroutines when computing 4p PSFs. Of course, this presumes subroutines which work irrespective of how their input operators have been generated and whether they are sliced or not. There is then no need to write separate routines for all those diagrams above depicting n>n or n>n refinements; instead, their contributions are generated automatically through nested function calls.

VI. FROM DISCRETE PARTIAL SPECTRAL FUNCTIONS TO CONNECTED CORRELATORS
The PSF computations described in the preceding sections lead to numerical results with the following structure. For a given permutation p, the PSF S p (ε) defined in Eq. (11) is represented as where ε and ε are ( −1)-tuples, the former containing the continuous variables ε i = ω 1···i (i < ), the latter consisting of discrete frequencies defined on the logarithmic grid used for binning and slicing, and The coefficients S p (ε) on the right, distinguished from the function S p (ε) on the left by having arguments with underbars, are stored as an ( −1)-dimensional histogram.
Such PSFs have to be convolved with a suitable kernel to compute p correlators. If the spectral representations of Sec. II are expressed through ε i variables, the convolution integrals d −1 ω p over functions of ω p (subject to ω 1··· = 0) turn into d −1 ε integrals over functions of ε, yielding ε sums upon insertion of Eq. (72) for S p (ε). For ZF correlators, e.g., we obtain [cf. Eqs. (3) and (4)] Similar expressions are obtained for MF and KF correlators. Evidently, each grid point ε defines the position of a pole for the kernel. To avoid a proliferation of poles, we choose the same ε grid for all permutations p.
Given a set of PSF histograms, the numerical evaluation of the above expressions is conceptually simple but requires great care in practice. One reason is that the ZF and KF kernels have poles lying very close to the real axis. The corresponding correlators, expressed as discrete ε sums over such kernels as in Eq. (73b), therefore do not yield smooth functions of frequency. To recover ZF or KF correlators with a smooth frequency dependence, as desired for a system whose original bath (before discretization) is continuous, the δ peaks of the PSFs have to be broadened. However, the broadening scheme should not spoil delicate cancellations in the sum p over permutations-a challenging requirement. By contrast, MF correlators need no broadening: their frequency dependence is smooth, since it enters as iω, so that the kernel denominators are always sufficiently far from zero. Vanishing bosonic frequencies in denominators are avoided, as their contributions are encoded in anomalous terms [cf. Eq. (8)].
The second major challenge arises if one is interested in computing the connected part, G con = G − G dis , of a correlator. It describes mutual correlations between particles (in contrast to their independent propagation), and is needed for extracting a corresponding vertex by amputating the external legs. Since the disconnected part G dis has very large extremal values, its subtraction has to be performed very accurately, and in such a manner that the discrete poles of G and G dis are aligned.
The present section describes how we deal with these challenges, with focus on computing G con for fermionic 4p correlators. We employ a multistep procedure: (A) Disconnected contributions are subtracted already at the level of PSFs. (B1) The PSFs are broadened in a Lorentzian manner and (B2) convolved with kernels to obtain G con ; in the MF, (B1) is omitted. (C) An equation-of-motion (EOM) scheme is used to improve the quality of G con ; this eliminates small disconnected contributions which remain after step (A) due to numerical inaccuracies. (D) External legs are amputated to obtain the vertex F . For the ZF and KF cases, the external legs are computed using the same Lorentzian broadening scheme as in (B1); improved results for G con can be obtained by subsequently reattaching external legs, computed using a more refined log-Gaussian broadening scheme. In the ensuing sections, we motivate and describe these steps in detail. Numerical results for G con are presented in Sec. VII below; results for F can be found in Ref. 1.

A. Subtracting disconnected parts
Any multipoint correlator can be split into a connected and a disconnected part, G = G con + G dis . For example, for ZF 4p correlators involving four fermionic operators, with each summand a product of two 2p correlators. The subtraction of G dis from the full G is numerically nontrivial for two reasons. First, the extremal values of G dis are generally much larger in magnitude than G con : In the ZF and KF, the 2p frequency conservation inherent in G dis is implemented via Dirac δ functions [cf. Eq. (74)] whose peak heights, even if regularized, are very large. In the MF, the Dirac δ's are replaced by Kronecker δ's, with a potentially large prefactor 1/T . Second, the δ's (Dirac or Kronecker) implementing 2p frequency conservation restrict the three-dimensional (3D) argument of G(ω) to a 2D plane. This 2D dependence is incompatible with the broadening scheme used in this work, which broadens 4p PSFs along all three frequency directions (see Sec. VI B below). Indeed, the 2D dependence is a peculiarity of the disconnected part; if a 4p correlator does not involve operators describing decoupled degrees of freedom, its connected part should have a full 3D structure. Given these challenges, it is advisable to subtract the disconnected part prior to broadening, already at the level of PSFs. To this end, we exploit the fact [1] that the decomposition G = G con + G dis has a counterpart at the level of PSFs, S p = S con Rather than computing the 2p PSFs occurring herein anew, we extract them from the (full) 4p PSFs via generalized sum rules, without extra NRG calculations. This facilitates the cancellation between S p and S dis p , as then all (4p and 2p) PSFs are obtained with the same accuracy. The sum rules are of the type These rules, and the (anti)commutation relation d σ d † σ − ζd † σ d σ = 1, imply that a 2p PSF can be obtained from a double integral over two 4p PSFs, e.g., As consistency checks on our numerical data, we verified that the 2p PSFs σ , d σ ](ε 1 ) = 1 with machine precision and (ii) match those obtained from direct 2p computations with an absolute accuracy of 10 −3 (see Fig. 14 in Appendix B). Property (i) is expected, since our iterative scheme for computing p PSFs conserves sum rules by construction. The quality of agreement for (ii) is very satisfying, given the complexity of the computation. Nevertheless, small discrepancies remain, as discussed in Appendix B. Therefore, it is preferable to compute S dis p using 2p PSFs obtained not directly, but from precisely the 4p PSFs involved in S p . Then, numerical inaccuracies inherent in S p are passed on to S dis p , facilitating the proper removal of disconnected (75), and using sums on ε 1 and ε 3 for the first and second factors in each line, we obtain We seek to express S con p = S p − S dis p in the form (72), with histograms S p (ε) and S dis p (ε) defined on matching 3D ε grids, so that their difference directly yields the histogram S con p (ε). (Different grids for S p and S dis p would yield nonmatching pole positions for G and G dis and a less reliable subtraction.) In Eq. (78), the first term has support in the ε 2 = 0 plane of the 3D grid. The second and third terms, however, have support at ε 2 = ε 1 + ε 3 , defining slanted planes incommensurate with the 3D grid. To associate these terms with points on the 3D grid, we rebin ε 1 + ε 3 onto the ε 2 grid, e.g., assigning a histogram element associated with the transition energies (ε 1 , ε 1 +ε 3 , ε i ) (i = 3, 1) to the grid point (ε 1 , ε 2 , ε i ) closest to it on a logarithmic scale (i.e., with ln |ε 2 | closest to ln |ε 1 + ε 3 |).
Since rebinning, needed to align the poles of G and G dis , slightly shifts some transition energies without recomputing the corresponding matrix elements, it introduces a non-negligible numerical error. Furthermore, at finite temperature, the effective finite length of the Wilson chain [22,71] introduces a low-energy cutoff, which limits the spectral resolution below T . More specifically, the thermal density matrix n DD is mostly populated at shells n ∼ n T = −2 log Λ (T /D) [22,71]. Therefore, transitions originating from very low-lying energy levels with E 1 T , resolved only within later shells with n > n T , are not captured accurately. Because of these numerical artifacts, arising from rebinning and finite temperature, the disconnected part in S is not canceled exactly. For example, Figs. 2(b) and 3(b) below show results for S con for a noninteracting system; it should vanish by Wick's theorem but does not. We remedy this by removing the remaining disconnected contributions with a EOM-based strategy, described in Sec. VI C below. Yet, before that, a broadening procedure, described next, is required for ZF and KF kernels.

B. Broadening
For ZF and KF correlators, whose kernels have poles lying very close to the real axis, the PSFs must be broadened before convolving them with kernels, as mentioned above. To this end, each Dirac δ(ε − ε) in S p (ε) is replaced by a peak-shaped broadening kernel, This amounts to replacing the convolution kernels, e.g., K(ω p , ε) in Eq. (73b), by broadened versions, For each grid direction i, the broadening kernel δ b (ε i , ε i ) is characterized by the shape and width of its peak when viewed as a function of ε i for given ε i . We first discuss the requirements for its width, say γ b,i . (The meaning of "width" of course depends on the peak shape, but the ensuing arguments are generic.) Because of the logarithmic bath discretization, the dominant contributions to the coefficients S p (ε) lie within clusters on the ε grid, with each cluster involving transitions within a specific combination of NRG shells. When plotted using logarithmic grid axes, the clusters are spaced by ∼ ln Λ, along each grid axis (see, e.g., Fig. 2 below), reflecting the single-particle level spacing of the discretized bath. The width of the broadening peak, plotted on the same axes, should be comparable to this in order to smooth out the clusters, i.e., ln (|ε i | + γ b,i )/|ε i | ln Λ. This requires a peak width proportional to frequency, γ b,i = b|ε i |, with b Λ − 1. Smaller choices for b are possible if one averages over several, say n z , discretization grids ("z averaging" [20]). We next discuss the choice of peak shape for δ b (ε, ε i ). NRG calculations of 2p functions conventionally use a log-Gaussian broadening kernel [20,22], i.e., a Gaussian function of ln |ε i /ε i | (for details, see Appendix C). It has the useful property of not overbroadening low-energy features resolvable only on a logarithmic scale.
For 4p correlators, we were not able to obtain satisfactory results using log-Gaussian broadening. By contrast, a Lorentzian broadening kernel with width b L |ε i |, yielded very good results for benchmark checks against various analytic results for the 4p vertex [1]. Via Eq. (79), it leads to broadened ZF and KF kernels of the form which match Eqs. (4) and (10), respectively, except for the replacements, (81c) Note that ξ 1···i and ξ j··· are the broadened version of ξ 1···i and ξ j··· , respectively, and not the sum of individual ξ m 's. Thus, Lorentzian broadening modifies the kernels only by shifting their poles away from the real axis by ±ib L |ε i |. Importantly, however, it does not change their analytic structure. As a result, the intricate cancellation patterns between different PSFs in permutation sums p largely remain intact during broadening. We believe this to be the main reason why, for 4p correlators, Lorentzian broadening outperforms log-Gaussian broadening.
The disadvantage of Lorentzian broadening, of course, is that Lorentzian peaks have long tails which very strongly overbroaden low-energy features-indeed, this is the reason why Lorentzian broadening is not used for 2p correlators. For 4p correlators, the effects of overbroadening appear to be more tolerable. Moreover, when computing the 4p vertex, some of their consequences can be canceled out by also using Lorentzian broadening for the externalleg 2p functions being divided out. To calculate the latter, we use the kernel (81) for = 2, but there is one subtlety: While we use the same value of γ i = γ 0 for all 1 ≤ ≤ 4 when computing 4p correlators, the values of γ i=1,2 for each 2p correlator depend on which leg it is attached to. This is explained in more detail in Sec. VI D.
The values used for the broadening parameters b L and γ 0 depend on the discretization parameter Λ, and on n z if z averaging is used. Artifacts appearing for frequencies lower than temperature can be smeared out by choosing γ 0 to be ∼ T . In practice, we choose b L and γ 0 just large enough to yield smooth results for broadened correlators. For the ZF and KF results in this paper, we used Λ = 4, n z = 4, γ 0 = 3T and b L = 0.6. (These values were also used for the results presented in Ref. 1, with two exceptions: we used b L = 1.4, γ 0 = T /2 for the x-ray-edge singularity, and b L = 0, γ 0 = T for the Hubbard atom. Though we treated the former in the ZF, we chose a small but finite T = 10 −5 D to define the density matrix and γ 0 .)

C. Equation of motion
When computing the connected contributions of PSFs, as described in Sec. VI A, some disconnected parts remain, due to numerical inaccuracies. To remove these, we employ an EOM strategy [74] which expresses G con through two 4p correlators whose disconnected parts mutually cancel. (This step is not needed if the computation of PSFs is exact, as was the case for benchmark results on the Hubbard atom in Ref. 1.) We first explain the EOM strategy in the MF, then discuss its transcription to the ZF and KF.
In NRG computations of 2p correlators for quantum impurity models, it is common practice to improve the results by computing the self-energy via an EOM scheme [75]. In the MF, QMC computations of local 4p correlators use a similar strategy, employing EOMs to obtain improved results for the connected part, G con [74,76]. Transcribed to our situation, Eq. (26) of Ref. 74 reads and d † σ are composite operators, and we put tildes on correlators containing one of them. Suppressing arguments to reveal the structure, this reads Equation (83) correctly yields G con 4p = 0 if H int = 0, since then G 4p = G 2p = 0 by construction. In the QMC study of Ref. 74, computing G con 4p via Eq. (83), rather than G 4p − G dis 4p , indeed led to markedly improved results. The reason is that the statistical errors for G 4p and G dis 4p are very different, whereas those for G 4p and G 4p are similar. Our NRG computations of G 4p and G dis 4p are also not perfectly matched-before subtracting the latter, rebinning of discrete spectral data is needed, see Sec. VI A. Thus, computing G con 4p via Eq. (83) is advisable in our case, too. Benchmark checks against analytical results (e.g., for weak interactions) show that a direct application of Eq. (83), with full 4p correlators as input on the right, does not yield optimal results for G con . Instead, we use connected 4p correlators as input, available via the subtraction scheme of Sec. VI A. Expressing both 4p correlators on the right of Eq. (83) as sums of connected and disconnected parts, we obtain as the disconnected parts mutually cancel, G 2p G dis 4p − G 2p G dis 4p = 0. This cancellation follows by using Eq. (74) and G d σ , d † σ = G d σ , d † σ . The latter equality stems from two different EOM derivations for G 2p , differentiating it with respect to either the first or the second time argument [74] and using ∂ t1 G 2p = −∂ t2 G 2p .
We evaluate Eq. (84) for G con 4p by inserting, on its right, results for G con 4p and G con 4p obtained by subtracting disconnected parts as described in Sec. VI A. Benchmark checks show that the output G con 4p is improved relative to the input G con 4p . The reason is likely that remnant disconnected contributions, still contained in the input correlators due to numerical inaccuracies, tend to cancel in Eq. (84).
We now turn to real-frequency EOMs for ZF and KF correlators. Apart from the fact that, here, input correlators must first be broadened, one proceeds analogously. The ZF, just as the MF, involves time-ordered correlators. Hence, ZF EOMs can be obtained from MF ones by simply replacing real by imaginary times. In the frequency domain, this yields the EOM (82) with iω replaced by ω.
In the KF in the Keldysh basis, a derivation analogous to that of Ref. 74, but tracking Keldysh indices, yields the same operator and frequency structure as Eq. (83): The only complication is that the ω 1 leg of each 4p correlator is contracted with the Pauli matrix σ . As a consistency check, we note that the ZF version of Eq. (83) can be obtained from Eq. (85) by transforming the latter to the contour basis and considering only components restricted to the forward branch of the Keldysh contour.
The mutual cancellation of disconnected parts discussed above applies to the real-frequency EOMs, too. For the ZF, Eq. (84) applies unchanged. For the KF, it becomes while the cancellation of disconnected parts can be verified The latter relation follows by evoking the fluctuation-dissipation theorem, G 22 (ω) = tanh(βω/2)(G 21 − G 12 )(ω), for both G 2p and G 2p .
The EOM scheme described above treats the frequency arguments of 4p correlators asymmetrically, in that ω 1 receives special treatment. For KF correlators G k , this leads to the fact that the k = 2111 component is obtained with better quality than others. This problem can be remedied using a symmetric EOM scheme presented recently in the MF in Ref. 76. We leave its transcription to the KF and use for NRG computations to the future.

D. Amputating and reattaching external legs
The strategies described in the preceding sections yield a connected 4p correlator, G con , free from disconnected parts. To obtain the corresponding 4p vertex F , the external legs are amputated by dividing out 2p correlators for all frequencies ω i . The MF vertex, e.g., is with a positive (negative) sign in the propagator of an outgoing (incoming) leg, associated with an annihilation (creation) operator in the definition of the 4p correlator.
The ZF and KF vertices are obtained using similar expressions (see Ref. 1 for details), but require special care. First, the choice of imaginary shifts for the arguments of the external legs G(ω i ) have to be chosen consistent with those of G con (ω) (see Appendix D of Ref. 1). We choose γ 1 = γ 0 and γ 2 = 3γ 0 (γ 1 = 3γ 0 and γ 2 = γ 0 ) for the outgoing (incoming) legs, where γ 0 is the value of all γ i 's in the kernel (81) used for computing the 4p correlator. Second, since real-frequency versions of G con are computed using Lorentzian-broadened kernels, we compute the 2p correlators in the denominator using the same Lorentzian broadening scheme (81c). Hence, the broadening width parameter b L is chosen identical to the one of the 4p correlator. This ameliorates the undesired overbroadening effects of Lorentzian broadening, in that they tend to cancel out in the ratio (87). This strategy is essential for obtaining the correct large-frequency behavior for F , since the large-frequency behavior for G con is dominated by the 2p correlators.
Having computed F , a yet-again improved version of G con can be obtained by multiplying F with externalleg 2p correlators, now computed using the customary log-Gaussian broadening scheme (cf. Appendix C). This ensures that those features dominated by the external legs are not overbroadened. In other words, to obtain a ZF or KF G con most accurately, we first compute the F vertex using Lorentzian-broadened ingredients, and then reattach the external legs through optimally broadened 2p correlators. This strategy for optimizing the treatment of external legs is particularly useful from the perspective that experimental probes typically measure response functions corresponding to correlators with external legs.

VII. RESULTS: CONNECTED CORRELATORS
To establish the power of our multipoint NRG scheme, we have performed detailed benchmark computations for the paradigmatic AIM. We presented MF and KF results for its 4p vertex in Ref. 1. To complement these, we here analyze both the underlying PSFs and the connected 4p correlators G con obtained from them by convolution with suitable MF, ZF, or KF kernels.

A. Model
The Hamiltonian of the AIM in standard notation is It describes a band of spinful electrons, hybridizing with a local level with energy d and Coulomb repulsion U .
We take d = −U/2. The coupling between impurity and bath is fully characterized by the hybridization function ∆(ν) = π|V | 2 δ(ν − ), which we choose box-shaped, ∆(ν) = ∆ θ(D−|ν|), with half-bandwidth D. Our goal is to compute the connected part of the local 4p correlator . For the NRG calculation, the bath is represented on a discrete, logarithmic grid, with grid intervals bounded by the points ±DΛ −k−z . Here, k ≥ 0 is an integer and z ∈ (0, 1] a shift parameter. A discrete energy representing each grid interval is chosen using the prescription proposed in Ref. 23. We use Λ = 4, and average results from n z different discretization grids, with z = 1/n z , 2/n z ,..., 1. We choose n z = 2 for MF correlators, which are less sensitive to discretization artifacts (for the same reason as why MF PSFs need not be broadened; see Sec. VI B), and n z = 4 for the more challenging real-frequency correlators. The Wilson chain is diagonalized iteratively by keeping 300 multiplets respecting U(1) charge and SU(2) spin symmetries. We generate and manipulate non-Abelian symmetric tensors using the QSpace library developed by Weichselbaum [77,78]. By exploiting non-Abelian symmetries, the different spin components of the 4p correlator are all obtained simultaneously.

B. Discrete partial spectral functions
The first output of our multipoint NRG calculations are discrete 4p PSFs. Figures 2(a)  and shown as a function of ε 1 and ε 3 at ε 2 = 0. Generally, we observe that the spectral contributions are spread over a range of energies, from −T to the largest energy scale in the system (here D). As 4p PSFs describe excitation spectra-their arguments are transition energies, ε = (E 21 , E 31 , E 41 )the spectral weight predominantly lies at positive energies. The spectral weight at negative energies, ε −T , is exponentially suppressed by the Boltzmann factor.
For noninteracting systems, Wick's theorem implies that 4p PSFs consist of only disconnected contributions, so that exact numerics would yield S = S dis , with (75). This factorization into two 2p PSFs is visible in Fig. 2(a), as the spectral weight is distributed on an equidistant, square grid for ln ε 1 and ln ε 3 , with grid spacing ln Λ. Though S con = S −S dis should vanish for U = 0, we numerically obtain a small but nonzero result, see Fig. 2(b). One reason stems from the fact that discrete spectral data at low energies, |ε i | T , are inaccurate. We therefore broaden them, using a broadening width ∼ T [24], to obtain S con 0 at energies |ε i | T . The finite but minuscule value of S con for |ε i | slightly larger than T [note the logarithmic color scale in Fig. 2(b)] comes from the (ε) (out of the 24 PSFs Sp), as a function of ε1 and ε3 at ε2 = 0. The top and bottom rows show noninteracting (U = 0) and strongly interacting (U/∆ = 5) cases, the left and right columns the full PSF S and its connected part S con , respectively. For U = 0, S con should vanish by Wick's theorem; the fact that small but nonzero values remain in (b) is due to numerical artifacts discussed in Sec. VI A. For U = 0, the distribution of spectral weight reflects the energy scales U/2 and TK D/200, and S con reveals mutual twoparticle correlations. These PSFs were computed using Λ = 4 and z = 1. For visualization, the discrete data was minimally broadened, using the log-Gaussian-Fermi kernel (C1) [24] with narrow width parameters. The intensity values for each row are normalized by the maximum magnitude of S on the left. tail of the inaccurate low-energy region, which is not fully smeared out through broadening.
Additional numerical artifacts arise from the rebinning of the sum of frequencies ε 1 + ε 3 into ε 2 bins during the computation of S dis (cf. Sec. VI A). This is exemplified in Fig. 3 (ε) as a function of ε 1 and ε 2 at ε 3 = ε 1 , with a layout analogous to Fig. 2. The rebinning artifacts, seen in Fig. 3(b) along the diagonal, are small as well. For example, the maximum magnitude of the data of Fig. 3(b) is about 2.5% of that of Fig. 3(a). Meanwhile, Fig. 2 is free from rebinning artifacts, because ε 2 = 0 there. Generally, rebinning artifacts are smeared out by broadening (cf. Sec. VI B). They can be further reduced via the EOM strategy described in Sec. VI C. (For U = 0, the property G con = 0 is then obtained exactly.) In the presence of interactions, the 4p PSFs exhibit a rich structure; see Figs. 2(c), 2(d), 3(c), and 3(d). The spectral distribution patterns in Figs. 2(c) and 3(c) change noticeably across characteristic energy scales induced by interactions, such as the Hubbard-band position U/2 and the Kondo temperature T K . The connected part S con , describing mutual two-particle correlations, now has con-L FIG. 3. Analogous to Fig. 2, but for a same-spin 4p PSF shown as a function of ε1 and ε2 at ε3 = ε1. In (b), the narrow adjacent red and blue regions along the diagonal ε1 ε2 > 0 reflect contributions from S and S dis which fail to cancel fully, due to rebinning involved for S dis . siderable weight over a wide frequency range, see Figs. 2(d) and 3(d). It is most pronounced at energies below U/2 and around T K but very small at large energies ε i > U ; hence, two-particle correlations are negligible there. The full and connected PSFs for the other permutations exhibit similar behavior; see Fig. 15 in Appendix D.

C. Connected 4p correlators for AIM
We next present results for the connected correlator G con , computed from PSFs using the strategies described in Sec. VI, including the amputation and reattachment of external legs (Sec. VI D). We start with MF results, as these can also be obtained by means of ED or QMC, and consider the AIM at strong interaction U/∆ = 5, with U/D = 1/5. In Ref. 1, we benchmarked the NRG 4p vertex, F σσ (iω), against QMC at an intermediate temperature T = 10 −2 D, finding relative deviations on the level of 1%. (For that benchmark, both QMC and NRG results were generated without using the EOM strategy of Sec. VI C.) We also presented results for the 4p vertex at much lower temperatures, out of reach for QMC algorithms. For convenience, Fig. 4(a) reproduces the results from Fig. 8 of Ref. 1 for F σσ (iω), now obtained by using the EOM strategy. Figure 4(b) shows the corresponding connected correlator G con σσ (iω) at T = 10 −4 D as a function of ν and ν at ω = 0, in the particle-hole representation [1,2], obtained by reattaching external legs to the vertex. These external legs ensure that G con decays in all directions (note the logarithmic color scale), contrary to the finite background value of F . There is a strong signal at the origin and along the axes ν ( ) = 0. Just as for the vertex [1], the equal-spin component vanishes for ν = ν and ω = 0 by symmetry. These features are mostly known already; Fig. 4 demonstrates that NRG provides a viable tool to compute MF 4p correlators at very low temperatures. For the real-frequency frameworks, ZF and KF, we cannot resort as easily to previous results. Hence, we next benchmark our NRG ZF and KF data against perturbation theory (PT), for the AIM at weak interaction, U/∆ = 1/2. In Ref. 1, we discussed second-order PT for the vertex in the MF and KF (one proceeds analogously in the ZF); results for G con follow after reattaching external legs. Figure 5(a) shows G con σσ (ω) for the ZF, again at T = 10 −3 D and ω = 0. There is very good agreement between NRG and PT. The color plots, now involving real and imaginary parts, show features somewhat similar to the MF results discussed above: a decay in all directions and a plus-shaped structure around a prominent signal at the origin. The two spin components mostly have opposite signs; however, contrary to the MF, the ZF data also involve sign changes within a spin component. Further, for these weak-interaction results, the ↑↓ component is larger in magnitude and extends to larger ν ( ) values than the ↑↑ component, since the former has a first-order contribution while the latter starts at order U 2 .
Turning to KF results, our scheme of extracting G con out of the full PSFs S works best for the retarded Keldysh components, where the Keldysh indices k of a correlator contain only one entry equal to 2. The reason is that we use a simple but asymmetric EOM (see Sec. VI C). Figure 5(b) shows the retarded correlator G k;con σσ (ω) with k = 2111 for the weakly interacting AIM (identical parameters as above). Again, there is good agreement between NRG and PT. Similar to the ZF results, we observe that the ↑↓ component exceeds the ↑↑ component in magnitude.
However, the structure of these ν-ν color plots at ω = 0 is markedly different between the ZF and KF: whereas the ZF ↑↓ component is mirror and point symmetric, the retarded counterpart is only point symmetric. The mirror symmetry of the ZF ↑↑ component is merely broken by the vanishing diagonal; the retarded counterpart, sharing the vanishing diagonal, is generally only point symmetric.
The reduced symmetry of the KF correlator is due to the fact the Keldysh indices k = 2111 single out the first index, related to the frequency ν. The ν-ν structure at ω = 0 of other retarded components not shown, i.e., k ∈ {1211, 1121, 1112}, follows from the k = 2111 component by 90 • rotations in either the real part, the imaginary part, or both. Finally, we turn to real-frequency connected correlators for the AIM at strong interaction, setting U/∆ = 5 as for the MF results. Figures 6(a) and 6(b) show the ZF and retarded KF (k = 2111) connected correlators G con σσ , respectively. We see that the symmetry considerations from above still hold, and the characteristic structure with distinct values at the origin, along the ν ( ) axes, and at the diagonal (with zero values for σ = σ ) persists. Apart from these similarities, the real-frequency structure seen in Fig. 6 is much richer than that of the weak-interaction results of Fig. 5 or the imaginary-frequency MF results of Fig. 4. The latter were obtained not only with identical system parameters, but from the same discrete data for the connected PSFs as the ZF and KF results of Fig. 6. This striking increase in complexity of real-frequency versus imaginary-frequency data illustrates the benefits of having direct access to real-frequency data-attempting to recover it from imaginary-frequency data via analytic continuation would be extremely challenging.
Altogether, the collection of plots in this section illustrates the power of the spectral representation of multi- In both panels, the real-frequency structure is much richer than that obtained at weak interaction (Fig. 5), and also than that of the imaginary-frequency MF results (Fig. 4) obtained for the same system parameters.

VIII. RESULTS: RIXS
We conclude our presentation of numerical results with an application of great physical interest: resonant inelastic x-ray scattering (RIXS). This type of measurement involves photon-in-photon-out spectroscopy. An incident x-ray photon excites an electron from a core level into a valence orbital [ Fig. 7(a)]. Then, a valence electron relaxes to fill the core hole, emitting a photon, while the difference between the energies of the incident and emitted photons, ω loss = ω in − ω out , is transferred into the solid [ Fig. 7(b)]. Measuring ω loss thus yields spectroscopic information about the excitations of the solid.
Below, we compute RIXS spectra by convolving a single 4p PSF for a specified set of four operators with a suitably chosen convolution kernel (different from those used hitherto). By using NRG, we achieve a fine resolution of power-law singularities and characteristic low-energy scales not accessible by ED methods.

A. Models
We have performed proof-of-principle computations for two minimal models. The first, to be called Mahan impurity model (MIM), was used by Mahan for describing x-ray absorption spectroscopy (XAS) in metals, manifesting the celebrated x-ray-edge singularity [70,[79][80][81] in the absorption spectrum. Its Hamiltonian is where c = c . The first term describes a flat, half-filled band of spinless electrons with ∈ [−D, D] and Fermi energy at = 0, the second a localized core level with energy p −D, usually occupied. An x-ray absorption process, described by the transition operator T † = c † p, transfers an electron from the core level into the conduction band. The resulting core hole, with hole number operator pp † , induces an attractive, local scattering potential −U p < 0 for the band (with U p | p |), described by the third term. Its sudden switch-on changes the wave functions of all conduction electrons, so that their initial and final ground states become orthogonal. This is Anderson's orthogonality catastrophe [82][83][84]. It is also responsible [85,86] for the singular behavior in the XAS spectrum [70,[79][80][81], accessible by NRG [86][87][88], and also to the RIXS spectrum [89,90]. To our best knowledge, the singularity-related features of the latter have not yet been studied numerically.
The second model describes x-ray absorption for an augmented Anderson impurity model (AAIM). We take the AIM Hamiltonian (88), augmented by core level terms, We ignore the spin of the core electron, since the interaction between the core and impurity levels is spin independent and the core-hole number cannot be larger than one. The transition operator describing x-ray absorption now reads T † σ = d † σ p. The XAS and RIXS spectra will have a richer structure than those for the first model, reflecting properties of the local density of states of the AIM, such as the presence of a Kondo resonance and Hubbard side bands.
We note in passing that a model of this type has been used to describe optical absorption for quantum dots with Kondo correlations [88,91], predicting x-ray-edgetype lineshapes that have been verified experimentally [92]. Moreover, our methodology is applicable also if the local density of states is computed self-consistently via DMFT, as done for RIXS computations for material systems [65,66]. We leave such applications, and generalizations to nonlocal RIXS processes with momentumdependent spectra, for future work. Here, our goal is to use the above two minimal models to highlight features in RIXS spectra that emerge when truly low-energy scales, beyond the reach of ED [55][56][57][58][59][60][61][62][63][64][65][66] or Bethe-Salpeter treatments [67][68][69], can be resolved.

B. Partial spectral functions for XAS and RIXS
We now turn to the computation of XAS and RIXS spectra, starting with the former. XAS measures the probability for the absorption of an x-ray photon [ Fig. 7(a)], the first part of a RIXS process, without tracking what happens subsequently. This probability is given by [54,93] where ω in + ω th is the incident photon energy, ω th the absorption threshold, and Γ p the inverse lifetime of the core hole. Using notation encompassing both models, we write the transition operator as T † σ = x † σ p, with x † = c † and the spin index omitted for the MIM, or for the AAIM. The transition operator links two Hilbert subspaces, to be called "no hole" and "one hole." These are not coupled by the Hamiltonian, [p † p, H] = 0, and can be diagonalized separately. The sums 2 or 1 run over all energy eigenstates of the one-hole or no-hole subspaces, respectively. The difference between their ground-state energies defines ω th . Equation (92) has the form of a convolution of a 2p PSF, S XAS σ (ε) = S[T σ , T † σ ](ε), with a Lorentzian peak: The x-ray absorption spectrum thus measures the 2p PSF S XAS σ (ω in + ω th ), broadened by Γ p . Next, consider RIXS spectra. They measure the probability for the absorption and subsequent reemission of x-ray photons with energies ω in + ω th and ω out + ω th , respectively, with ω loss = ω in −ω out (Fig. 7). This probability is given by the Kramers-Heisenberg formula [53,54,94], In the second line, the expression within the square is the amplitude for a second-order process, involving upand downward transitions from the no-hole to one-hole subspaces and back. The third line explicates the dependence on the frequencies ω in and ω out . Transitions with ω in = E 41 − ω th and ω out = E 23 − ω th conserve energy and are "real." If both equalities hold (at T = 0 this requires ω in > 0), the net result is a resonant fluorescencetype process with a diverging amplitude due to vanishing energy denominators, requiring regularization by Γ p . If one or both equalities are violated, the net result is an off-resonance Raman-type process involving virtual intermediate states. Nevertheless, its amplitude can be sizable provided that ω loss = E 31 (at T = 0 this requires ω loss > 0), corresponding to a real initial-to-final transition between two no-hole states. Hence, I RIXS has support for both positive and negative ω in , but only for ω loss −T , i.e., the dependence on ω loss has threshold character. Sharp features in the ω in dependence of I RIXS are smeared out by Γ p . By contrast, the ω loss dependence enters separately from Γ p and ω in in Eq. (94); hence, it is not smeared out and can contain sharp, distinctive features even if ω loss Γ p , |ω in |. We can express Eq. (94) as a convolution involving a 4p PSF, S RIXS σσ (ε) = S[T σ ,T † σ ,T σ ,T † σ ](ε) and suitable kernel: .
The PSF arguments ε 1 and ε 3 are connected to ω in + ω th by Lorentzian kernels, ε 2 to ω loss by a Dirac δ function.
The dependence of I RIXS on ω in and ω loss thus probes that of S RIXS on ε 1 , ε 3 , and ε 2 , respectively, the former two broadened by Γ p , the latter not. [Note that, if Γ p → 0 + , I RIXS can also be viewed as a component of the KF correlator in the contour basis, (ω), with c = ++−− and frequency arguments, in the particle-hole convention (89), given by ν = ν = ω in +ω th , ω = −ω loss .] For sufficiently large |ω in |, the dependence of I RIXS on ω loss mimics that of a dynamical susceptibility. This can be seen as follows. When |ω in | is much larger than Γ p and also lies far outside the support of the PSF along the ε 1 − ω th and ε 3 − ω th axes, the denominator of the RIXS kernel in Eq. (95) can be approximated by |ω in | 2 , ignoring its ε 1 and ε 3 dependence. According to the PSF sum rules (76), I RIXS then reduces to a 2p PSF, The second line, following via T σ T † σ = x σ x † σ p † p, involves a 2p PSF evaluated purely within the no-hole subspace. It represents the imaginary part of an equilibrium 2p correlator, the ZF particle-hole susceptibility χ σσ (ω loss ) = G[x σ x † σ , x σ x † σ ](ω loss ), considered for frequencies ω loss T . [The latter restriction is needed because χ(ω) involves PSFs for two permutations, whereas Eq. (96b) involves only one.] Similarly, when Γ p , instead of |ω in |, is the largest scale in the kernel denominator, we have We next discuss some computational details. As mentioned above, the core-hole number is conserved; hence, the no-hole and one-hole subspaces can be diagonalized separately [71,84,86]. This yields two sets of AS bases, capable of accurately resolving low frequencies |ω in |, |ω loss | D. Their direct sum yields an eigenbasis for the whole Hilbert space, suitable for computing PSFs as explained in Sec. V. The density matrix is populated only by no-hole states. We write the eigenenergies of the one-hole subspace as E i + ω th . This splits off ω th from the arguments ε of S XAS (ε) and ε 1 and ε 3 of S RIXS (ε), canceling the ω th in the kernels of Eqs. (93) and (95). The discrete PSFs S XAS (ε) and S RIXS (ε) are saved as histograms. To focus on the inelastic part of I RIXS (ω loss 0), we exclude the spectral weights of the histogram S RIXS (ε) at ε 2 = 0. We use Λ = 4 and n z = 4. Broadened PSFs S XAS (ε) and S RIXS (ε) are obtained using centered log-Gaussian kernels, as described at the end of Appendix C. We have verified that our NRG code and a recent ED code [66] give consistent result for Wilson chains short enough that exact diagonalization is possible. Figure 8 shows a typical XAS spectrum for the MIM. I XAS (ω in ) features a threshold at ω in = 0, and, just above it, the x-ray-edge singularity [ Fig. 8(a)]. Sharp features are smeared out by nonzero Γ p or T . As these have analogous broadening effects, we take Γ p > T throughout, except for the scaling analysis in the very next paragraph and Fig. 9. In the following, we focus on the behavior outside the broadening region, ω in Γ p . A log-log plot for positive frequencies [ Fig. 8(b), black line] reveals powerlaw behavior in the range Γ p ω in 0.5D. Its form is consistent with the analytical prediction [70,[79][80][81]

C. MIM: XAS and RIXS results
and previous NRG studies [86,88]. Here, δ is the phase shift for conduction electrons near the Fermi level induced by the core-hole scattering potential. It can also be computed using δ = π∆ h , where ∆ h is the charge drawn in toward the scattering site by the core hole [84,86]. The power-law behavior can be traced back to Anderson orthogonality [82,83], see Ref. 86 and Appendix E for a heuristic discussion. Such power-law behavior is inacces-sible to ED methods, since the description of Anderson orthogonality in essence requires an infinite bath. Figures 9(a) and 9(b) show color-scale plots of the RIXS spectrum, using (a) linear and (b) logarithmic frequency axes. I RIXS (ω in , ω loss ) has support only at positive ω loss , and at both positive and negative ω in , with somewhat more weight at positive than negative ω in . It shows a power-law dependence on both ω in and ω loss , culminating in a divergence around ω in = 0 (cut off by Γ p ; to better reveal it, we here chose Γ p = 0 + ). The powerlaw exponents depend on the relative values of |ω in | and ω loss , as indicated by legends in Fig. 9(b). They were identified by plotting I RIXS as a function of ω loss or ω in , respectively, for several fixed values of the other variable [Figs. 9(c) and 9(d)]. We find (i) I RIXS ∼ |ω in | −2−2α ω loss for ω loss |ω in |; (ii) I RIXS ∼ |ω in | −1−2α for ω loss |ω in |; and (iii) I RIXS ∼ |ω in | −1−α ω −α loss for ω loss |ω in |, with the same exponent α as for XAS spectra. In summary, The function f is not fully symmetric, f (x) = f (−x), having the same asymptotic behavior but different prefactors for arguments of opposite sign. Figure 9(c) also shows that I RIXS (ω in , ω loss ), viewed as a function of ω loss for fixed ω in , has a maximum around ω loss |ω in |. Its value at this maximum, I RIXS max (ω in ), scales as ∼ |ω in | −1−2α [ Fig. 8(b), blue lines], with a smaller prefactor (by about an order of magnitude) for negative than positive ω in .
The XAS and RIXS power laws for the MIM are governed by the same exponent α because both stem from Anderson orthogonality. This was recognized already in the early 1970s [89], leading to analytic predictions for RIXS exponents [cf. Eqs. (37) and (41) These predictions match our numerics for case (i) above (numerically most challenging since ω loss |ω in |), but not for cases (ii) and (iii). The discrepancies in the latter two cases, where ω loss is not |ω in |, suggest that, in Ref. 90, the treatment of excitations in the no-hole subspace [labeled 3 in Eq. (94)] was not sufficiently accurate. Clarifying this would require revisiting the analysis of Ref. 90, which is, however, beyond the scope of this work.
Because of the strong dependence of the RIXS spectrum on ω in , it is convenient, when analyzing its dependence on ω loss , to first normalize the spectrum by its maximal value, I RIXS norm (ω in , ω loss ) = I RIXS (ω in , ω loss )/I RIXS max (ω in ). A colorscale plot of I RIXS norm is shown in Fig. 10(a). It highlights the presence of a peak as a function of ω loss for each value of ω in . We already encountered this peak in Fig. 9(c). Its evolution with ω in is analyzed in Fig. 10(b), showing very similar, though slightly asymmetric behavior for positive and negative ω in . For |ω in | in the range between Γ p ω in 0.5D, I RIXS norm has a peak at ω loss |ω in |. For positive or negative ω in , this implies that emission is strongest for ω out 0 or ω out 2ω in , respectively. In the former case, all the above-threshold incident energy is absorbed by the Fermi sea, implying fluorescence-type behavior. The latter case does not seem to have a simple interpretation, but for brevity we will call it "fluorescencelike," too. For both cases, the observed behavior implies strongly energy-dependent matrix elements T ij in Eq. (94), reflecting Anderson orthogonality.
For incident energies smaller than the hole decay rate, |ω in | < Γ p , the peak gets pinned to ω loss Γ p [ Fig. 10(b)], since the dependence on ω in becomes smeared. For very large incident energies, ω in 0.5D, the peak positions become pinned at ω loss D. To understand the latter, we evoke single-particle arguments, valid at large energies: The probability distributions for the excitation of a single particle or hole are flat (following the local density of states), with support in the range [0, D]. The distribution of a one-particle-one-hole excitation energy has range [0, 2D], with a peak around D.
This case study of the RIXS spectrum of the MIM shows (i) that our scheme is capable of resolving twovariable power-law behavior even if the two variables differ by orders of magnitude and (ii) that the normalized RIXS spectrum is dominated by a peak at ω loss |ω in | (provided that |ω in | > Γ p ), reflecting strongly energydependent transition matrix elements.

D. AAIM: XAS and RIXS results
Finally, we discuss our XAS and RIXS results for the AAIM. The unnormalized XAS and RIXS spectra are shown in Figs. 11 and 12, respectively. For |ω in | T K and ω loss T K , they-remarkably-exhibit power-law behavior of the same functional form as for the MIM.
For the XAS spectrum, this is expected: Its power-law form can be understood [85,86] by assuming Fermi-liquid excitations in the low-energy regimes of the no-hole or one-hole subspaces, and this assumption is valid for both the MIM and the AAIM. (Low-energy means ε ≤ D for the MIM, and ε ≤ T K for the AAIM, where T K is the Kondo temperature in the no-hole subspace.) The only difference is that, for the AAIM, the power-law exponent governing the XAS spectrum depends on spin, where ∆n dσ = n one dσ −n no dσ is the difference in n dσ = d † σ d σ , the dσ-level occupancy, between the one-hole and no-hole ground states. Equation (100)    exponent α close to 0.5 was viewed as a fingerprint of Kondo correlations in absorption spectra.
While the XAS spectrum in Fig. 11 matches the findings of Refs. 91 and 92, the RIXS spectrum of the AAIM, shown in Fig. 12, is new. As mentioned above, it shows power laws with the same functional form as for the MIM, Fig. 9. Presumably, this commonality again reflects the Fermi-liquid nature of the low-energy excitations of the MIM and AAIM. Understanding this in detail is left as a challenge for future work. From a purely computational perspective, Fig. 12 illustrates that our approach, applied to a model with strong interactions and subtle correlations, is able to probe energy regimes differing by many orders of magnitude, uncovering power-law behavior with a quantitatively consistent value for the exponent α.
Though the unnormalized RIXS spectra of the MIM and AAIM look similar, their normalized forms exhibit differences. For the AAIM, we define I RIXS norm,σσ = I RIXS σσ /I RIXS max , with I RIXS max (ω in ) = max ω loss ,σσ [I RIXS σσ (ω in , ω loss )]. As an example, Fig. 13 shows I RIXS ↑↑ . Compared to Fig. 10 for the MIM, it has a richer structure, reflecting the presence of Kondo correlations. Consider first the case that the inverse core-hole lifetime is much smaller than the Kondo temperature, Γ p T K [Figs. 13(a) and 13(b)]. At low incident energies, Γ p < |ω in | < T K , the spectrum shows fluorescencelike behavior, having a single peak at ω loss |ω in |. If |ω in | is further increased above T K , a Raman-type peak remains pinned at T K , more pronounced for I RIXS norm,↑↓ than I RIXS norm,↑↑ . Moreover, for positive ω in (but not for negative ω in ) a second, shoulder-type structure emerges, moving upward with ω in , in fluorescencelike fashion. (See the curves for ω in /D = 10 −1 , 10 −0.5 .) The height of this shoulder decreases as ω in gets larger, and for ω in D, the shoulder essentially disappears.
Next, consider the case Γ p T K [Figs. 13(c) and 13(d)]. Then, a fluorescence-type peak is found only at positive incident energies, ω loss ω in , and, since it requires Γ p ω in D, only for ω in well above T K ; in fact, such a structure is visible only as a weak shoulder on the right of Fig. 13(d). Importantly, however, the Raman-type peak at T K persists-its position and shape are the same as for Γ p T K [ Fig. 13(b)]. The reason is that ω loss enters the Heisenberg-Kramers formula separately from Γ p , as mentioned earlier.
The pinning of the peak at T K for large |ω in | is consistent with the relation (96b) between I RIXS and dynamical susceptibilities. Expressed through d operators, it reads (ωin, ω loss ) for the AAIM (same parameters as for Fig. 11). The layout is the same as for Fig. 9. The spectrum I RIXS

↑↓
(not shown) is qualitatively similar, exhibiting the same power laws.
local spin fluctuations, which has a peak at T K [95]. By contrast, the same-spin spectrum I RIXS ↑↑ receives contributions from both the spin and charge susceptibilities, and hence has peaks at both T K and U/2. The peak height of the latter is very much smaller than that of the former, since the spin and charge susceptibilities have peak heights inversely proportional to the peak positions.
A more detailed discussion of the RIXS spectra of the AAIM will be published elsewhere. For present purposes, the main messages of the proof-of-principle data presented here are the following: (i) at large |ω in | or large Γ p , RIXS spectra probe dynamical susceptibilities; (ii) for models involving Kondo correlations, they exhibit a distinct Raman-type peak at ω loss T K , well-resolved irrespective of the relative size of Γ p and T K ; and (iii) our multipoint NRG scheme is very well suited for uncovering such peaks, even for widely separate scales, ω loss |ω in |.

A. Summary
In this work, we showed how NRG can be used to compute local 3p and 4p correlators of quantum impurity models. Building on the spectral representation introduced in the accompanying paper, Ref. 1, we used NRG to compute the fundamental PSFs, and then convolved these with suitable kernels to obtain both imaginary-and real-frequency correlators.
To compute the PSFs, we first developed a refined tensor-network diagrammatic notation and explained how to expand general operators along the NRG Wilson chain in the AS basis. Then, we introduced a "slicing" technique that enables recursively expressing higher-through lowerpoint PSFs. The resulting PSFs can be computed at any temperature and contain spectral information over a wide range of energies, from the bare energy scales (e.g., Coulomb repulsion, bandwidth) down to excitations even below emergent energy scales (such as the Kondo temperature).
Once the PSFs have been obtained, imaginary-and real-frequency correlators can be computed from the same PSFs by applying appropriate convolution kernels. This is advantageous compared to other methods that mainly treat imaginary-frequency correlators, since the numerical analytic continuation of multipoint objects from imaginary to real frequencies is extremely challenging. Like the PSFs, the correlators can be computed at arbitrary temperature and frequencies. Our framework also encompasses a suite of strategies, including real-frequency EOMs, to accurately obtain connected correlators and vertex functions. Altogether, results of the former, shown here, and of the latter, shown in Ref. 1, pass numerous qualitative and quantitative benchmark tests.
In addition to standard fermionic 4p functions, we also applied our method to calculate RIXS spectra. To this end, we rephrased the Kramers-Heisenberg formula, a theoretical description of RIXS spectra, as a special convolution of a PSF involving transition operators. By studying two minimal models, we could identify distinctive contributions to the RIXS spectra: (i) a fluorescencetype peak at ω loss |ω in |, seen if ω in > Γ p , or if −ω in is larger than Γ p and smaller than the upper bound to (quasi)particle excitation energies (half-bandwidth D for the bare particles in the MIM; Kondo temperature T K for TK in (c) and (d). For both cases, the same layout is used as for Fig. 10: (a),(c) Color-scale plots using linear frequency scales, and (b),(d) line cuts using a logarithmic frequency scale to show the ω loss dependence for several fixed values of ωin, indicated by dots placed at ω loss = |ωin|. A fluorescence-type peak or shoulder is seen at ω loss |ωin| if Γp |ωin| TK, and at ω loss ωin if max{TK, Γp} ωin D. Moreover, a Raman-type peak occurs at ω loss TK for max{|ωin|, Γp} TK. the quasiparticles in the AAIM); and (ii) a Raman-type part, peaked at ω loss T K for the AAIM, seen if |ω in | T K or if Γ p T K . This Raman-like part is not smeared out by large Γ p , demonstrating the suitability of RIXS for probing (potentially small) many-body correlation scales. Our approach, which can treat strong correlations and low-energy excitations accurately, may provide a benchmark for more approximate, but computationally cheaper, methods.

B. Outlook
Our scheme is both general and powerful, lending itself to various interesting applications. The local 4p vertex, which we analyze in detail in Ref. 1 and here use to improve the accuracy of the connected 4p correlators, is at the heart of Feynman-diagrammatic analyses of strongly correlated systems. It is an input for diagrammatic ex-tensions of DMFT aiming to incorporate nonlocal correlations [2]. Real-frequency applications along these lines are discussed in Sec. V B of Ref. 1. Furthermore, the vertex is an input to the computation of nonlocal response functions, including the magnetic structure factor [96] and transport properties [3]. For the latter, it was recently found that vertex corrections can be sizable even if the self-energy is practically local, i.e., consistent with the DMFT approximation [97]. Since a major complication for computing conductivities is the numerical analytic continuation, our real-frequency method promises to be a breakthrough tool in this regard.
Another interesting application is to scrutinize how RIXS spectra reflect emergent energy scales in correlated lattice systems. Such scales indicate key physical mechanisms: For instance, the Mott metal-to-insulator transition in the Hubbard model is accompanied with the decrease of the spin screening scale [25,98,99], and Hund metals are characterized by spin-orbital separation [32, 35-37, 100, 101], i.e., the fact that spin and orbital screening scales are separated by orders of magnitude. Such small, many-body energy scales are hard to capture with established RIXS methods such as ED, while our NRG scheme is ideally suited for them.
All applications will benefit from further methodological development. The asymmetric EOM used in this work leads to significant improvement for only some components of Keldysh correlators. For applications building on the vertex, improvement for all Keldysh components, through a KF analogue of the symmetric EOM scheme of Ref. 76, will be necessary. Moreover, refined multipoint broadening schemes, rather than the Lorentzian broadening prone to overbroadening used here, should be devised to better resolve all multiparticle spectral features. Using Eq. (A1), all operator products can be reduced to nested sums over shells, n n≥n n≥n ..., with each summand involving only the building blocks of Eq. (A1), but containing numerous shell-off-diagonal products. For the computation of thermal averages using Eq. (48) though, we prefer a single sum n , involving no shell-offdiagonal contributions. This can be achieved by systematically using Eq. (36b), n>n =KK XX O n XX = O n KK , etc., to collect all shell-off-diagonal contributions into KK sectors of earlier shells. The result is a sum n containing, for given n, a sum over all possible products of shell-n projections, excluding the all-K case. The latter is represented, in refined fashion, via later shells with larger n. This scheme readily reproduces Eqs. (50). For example, expanding AB using Eq. (37), simplifying via Eq. (A1), and collecting shell-off-diagonal terms using Eq. (36b), we obtain which equals Eq. (50a). In the second line, the first two terms contain the nonzero shell-diagonal n = n contributions; the third term collects the shell-off-diagonal contributions with n > n into A n KK , and the fourth those with n > n into B n KK . Equations (50b) and (50c) can be found similarly.
However, with the 4p PSFs obtained by NRG, there exist small but finite differences among S 12;σσ and S 34;σσ . In Fig. 14, we compare them with a 2p PSF S 2p obtained by using our 2p method [cf. Eq. (59)], which is equivalent to fdm-NRG. While S 34;σσ is identical to S 2p up to double precision 10 −16 regardless of σσ (hence not shown), S 12;σσ differs from S 2p , and also from each other depending on σσ . The differences are smaller than the local maxima for individual clusters of spectral weights (which appear as peaks in the plot) by at least two orders of magnitude.
There are two reasons for these differences. First, the reduced density matrices in the kept sectors, n KK , are not diagonal, since they are obtained by tracing out entangled states over a subsystem. Second, the way of measuring the first argument ε 1 of 4p PSFs depends on the sectors. For example, while the first operator A is sliced with respect to ε 1 for computing S n KDKK (ε) [cf. Eq. (71c)], ε 1 is measured as the energy difference across the product A (not just A) for DX X X and KDDD [cf. Eq. (69)]. By contrast, ε 3 is always associated with D, up to a sign. For 2p PSFs, the argument ε is always associated with A or equivalently, B, up to a sign. Given that the n KK 's are generally not diagonal, such a discrepancy in measuring ε 1 introduces a numerical artifact to the distribution of spectral weights along the ε 1 grid, for both the full 4p PSFs and S 12;σσ (ε = ε 1 ). For S 34;σσ (ε = ε 3 ), the ε 1 dependence has been summed over, so the artifact is removed.
As mentioned in Sec. VI A, the connected part can have smaller values than the full 4p correlator, and thus be more susceptible to numerical noise. Hence, we need to remove the above-mentioned numerical artifact when subtracting the disconnected part. For this, in evaluating Eq. (75), we substitute the 2p PSFs obtained by summing over (ε 2 , ε 3 ), such as S 12;σσ , to the 2p PSFs involving O 1 , and those obtained by summing over (ε 1 , ε 2 ), such as S 34;σσ , to the rest. Moreover, to calculate the disconnected part for a given spin configuration (σ = σ or σ = σ ), we use the 2p PSFs derived from 4p PSFs for that configuration.
Bosonic correlators also arise when computing the XAS and RIXS spectra of Sec. VIII of this paper. We broadened them as follows. Energy dependencies not affected by Γ p were broadened using δ CL+F , with b F = T /5. This applies to the dependence of the RIXS kernel on ε 2 regardless of Γ p . For Γ p < T , it also applies to the dependence of the XAS kernel on ε, and of the RIXS kernel on ε 1 and ε 3 , because then their Lorentzian dependence on Γ p becomes irrelevant. By contrast, for Γ p > T , we implemented Lorentzian broadening of the ε or ε 1 and ε 3 dependencies by replacing δ F in Eq. (C1) by a Lorentzian δ L of the form (80), with b L |ε i | there replaced by Γ p . For the RIXS kernel, this yields the imaginary parts of its Lorentzian denominators; their real parts were then obtained using Kramers-Kronig transformations. We used the broadening parameters b CL = ln Λ, except for the data shown in Fig. 13, where a smaller value, b CL = (1/2) ln Λ, was used to better resolve the separation of the Raman-and fluorescence-type peaks.