Spectral representation of Matsubara n-point functions: Exact kernel functions and applications

In the field of quantum many-body physics, the spectral (or Lehmann) representation simplifies the calculation of Matsubara n-point correlation functions if the eigensystem of a Hamiltonian is known. It is expressed via a universal kernel function and a system- and correlator-specific product of matrix elements. Here we provide the kernel functions in full generality, for arbitrary n, arbitrary combinations of bosonic or fermionic operators and an arbitrary number of anomalous terms. As an application, we consider bosonic 3- and 4-point correlation functions for the fermionic Hubbard atom and a free spin of length S, respectively.


Introduction
Multi-point correlation functions of n quantum mechanical operators, also known as n-point functions, are a central concept in the study of quantum many-body systems and field theory [1].They generalize the well-known 2-point functions, which, for the example of electrons in the solid state, are routinely measured by scanning tunneling spectroscopy or angleresolved photon emission spectroscopy [2].For magnetic systems, the 2-point spin correlators can be probed in a neutron scattering experiment.Higher order correlation functions with n = 3, 4, 5... can for example be measured in non-linear response settings [3].In the emerging field of cold atomic quantum simulation, (equal-time) n-point functions are even directly accessible [4].
On the theoretical side the study of higher order correlation functions gains traction as well.One motivation is the existence of exact relations between correlation functions of different order n [5,6].Although these exact relations can usually not be solved exactly, they form a valuable starting point for further methodological developments like the parquet approximation [7].Thus even if the 4-point correlator (or, in that context, its essential part, the one-line irreducible vertex [1]) might not be the primary quantity of interest in a calculation, it appears as a building block of the method.Another example is the functional renormalization group method (fRG) in a vertex expansion [8,9].It expresses the many body problem as a hierarchy of differential equations for the vertices that interpolate between a simple solvable starting point and the full physical theory [10].Whereas experiments measure correlation functions in real time (or frequency), in theory one often is concerned with the related but conceptually simpler versions depending on imaginary time [1].In the following, we will focus on these Matsubara correlation functions, which, nevertheless feature an intricate frequency dependence.
Whereas the above theoretical methods usually provide only an approximation for the npoint functions, an important task is to calculate these objects exactly.This should be possible for simple quantum many body systems.We consider systems simple if they are amenable to exact diagonalization (ED), i.e. feature a small enough Hilbert space, like few-site clusters of interacting quantum spins or fermions.Also impurity systems, where interactions only act locally, can be approximately diagonalized using the numerical renormalization group [11].
Knowing the exact n-point functions for simple systems is important for benchmark testing newly developed methods before deploying them to harder problems.Moreover, n-point functions for simple systems often serve as the starting point of further approximations like in the spin-fRG [12][13][14], or appear intrinsically in a method like in diagrammatic extensions of dynamical mean field theory [15] with its auxiliary impurity problems.Another pursuit enabled by the availability of exact n-point functions is to interpret the wealth of information encoded in these objects, in particular in their rich frequency structure.For example, Ref. [16] studied the fingerprints of local moment formation and Kondo screening in quantum impurity models.
In this work we complete the task to calculate exact n-point functions by generalizing the spectral (or Lehmann) representation [1,17] for Matsubara n-point correlation functions to arbitrary n.We assume that a set of eigenstates and -energies is given.Following pioneering work of Refs.[18][19][20] and in particular the recent approach by Kugler et al. [21], we split the problem of calculating imaginary frequency correlators into the computation of a universal kernel function and a system-and correlator-specific part (called partial spectral function in Ref. [21]).We provide the kernel functions in full generality for an arbitrary number n of bosonic or fermionic frequencies.Previously, these kernel functions were known exactly only up to the 3-point case [18], for the fermionic 4-point case [19][20][21] or for the general n-point case [21] but disregarding anomalous contributions to the sum that the kernel function con- sists of.These anomalous contributions are at the heart of the complexity of Matsubara n-point functions.They occur when certain combinations of eigenenergies and external frequencies vanish individually, see the anti-diagonal rays in Fig. 1(c).Physically, they correspond to longterm memory effects, are related to non-ergodicity and, in the case of bosonic two-point functions reflect the difference between static isothermal susceptibilities and the zero-frequency limit of the dynamical Kubo response function [22,23].The structure of the paper is as follows: In Sec. 2 we define the Matsubara n-point function G A 1 ...A n ω 1 , ..., ω n−1 and review some of its properties.The spectral representation is derived in Sec. 3 with Eq. ( 15) being the central equation written in terms of the kernel function K n (Ω 1 , ..., Ω n−1 ).Our main result is an exact closed-form expression of this most general kernel function which is given in Sec. 4. Examples for n = 2, 3, 4, 5 are given in Sec. 5 where we also discuss simplifications for the purely fermionic case.We continue with applications to two particular systems relevant in the field of condensed matter theory: In Sec. 6, we consider the Hubbard atom and the free spin of length S, for which we compute n-point functions not previously available in the literature.We conclude in Sec. 7.

Definition of Matsubara n-point function
We consider a set of n = 2, 3, 4, ... operators {A 1 , A 2 , ..., A n } defined on the Hilbert space of a quantum many-body Hamiltonian H.The operators can be fermionic, bosonic or a combination of both types, with the restriction that there is an even number of fermionic operators.As an example, where d † and d are canonical fermionic creation and annihilation operators.A subset of operators is called bosonic if they create a closed algebra under the commutation operation.They are called fermionic if the algebra is closed under anti-commutation, see Sec. 1 of Ref. [24].Spin operators are thus bosonic.
We define the imaginary time-ordered n-point correlation functions for imaginary times τ k ∈ [0, β], [25,26], where [21] multiplies with (−1) n−1 .In Eq. ( 1), the imaginary time-ordering operator T orders the string of Heisenberg operators, where p is the permutation p ∈ S n such that τ p(1) > τ p(2) > ... > τ p(n) [see Fig. 1 Imaginary time-ordered correlation functions (1) fulfill certain properties which we review in the following, see e.g.[26] for a more extensive discussion.First, they are invariant under translation of all time arguments, with τ ∈ such that τ k + τ ∈ [0, β].They also fulfill periodic or anti-periodic boundary conditions for the individual arguments τ k , where ζ k = +1 or −1 if A k is from the bosonic or fermionic subset of operators, respectively.This motivates the use of a Fourier transformation, where 2)/β with m k ∈ are bosonic or fermionic Matsubara frequencies, respectively, and ω k is shorthand for m k ∈ .Note that fermionic Matsubara frequencies are necessarily nonzero, a property that will become important later.As we will not discuss the real-frequency formalisms, we will not write the imaginary unit in front of Matsubara frequencies in the arguments of G A 1 ...A n (ω 1 , ..., ω n ).Again, note that in the literature, different conventions for the Fourier transformation of n-point functions are in use.In particular some authors pick different signs in the exponent of Eq. ( 7) for fermionic creation and annihilation operators, or chose these signs depending on operator positions.Time translational invariance (4) implies frequency conservation at the left hand side of Eq. (7), where on the right hand side we skipped the n-th frequency entry in the argument list of G.Note that we do not use a new symbol for the correlation function when we pull out the factor β and the Kronecker delta function.

Spectral representation of G A 1 ...A n ω 1 , ..., ω n−1
The integrals involved in the Fourier transformation ( 7) generate all n! different orderings of the time arguments τ k .As in Ref. [21] it is thus convenient to use a sum over all n! permutations p ∈ S n and employ a product of n−1 step-functions θ , with θ (x) = 1 for x > 0 and 0 otherwise, to filter out the unique ordering for which To expose explicitly the time dependence of the Heisenberg operators, we insert n times the basis of eigenstates and -energies of the many-body Hamiltonian H. Instead of the familiar notation j 1 , j 2 , ... and E j 1 , E j 2 , ... we employ 1 , 2 , ... and E 1 , E 2 ,.... for compressed notation and denote operator matrix elements as A 1 2 = 1|A|2 .We obtain and apply the Fourier transform according to the definition (7), p(1) τ p(1) where we defined In Eq. ( 11), the first line carries all the information of the system and the set of operators {A 1 , A 2 , ..., A n }.The remaining terms can be regarded as a universal kernel function defined for general {Ω 1 , Ω 2 , ..., Ω n } probed at Ω k ∈ which depends on the system and correlators via (12).Upon renaming the τ-integration variables τ p(k) → τ k , this kernel function is written as follows: In the second line we split K n into a part K n proportional to βδ 0,Ω 1 +Ω 2 +...+Ω n and the rest R n .We dropped Ω n from the argument list of K n which can be reconstructed from {Ω 1 , ..., Ω n−1 }.Finally, we express G A 1 ...A n ω 1 , ..., ω n of Eq. ( 11) using the kernel K n so that the general k of Eq. (12).For these, Ω since the E k cancel pairwise.The structure of Eq. (8) (which followed from time translational invariance) implies that the terms proportional to R n are guaranteed to cancel when summed over permutations p ∈ S n , so that only the terms proportional to K n remain.We drop the βδ 0,ω 1 +ω 2 +...+ω n from both sides [c.f.Eq. ( 8)] and find the spectral representation of the npoint correlation function in the Matsubara formalism, An equivalent expression was derived in the literature before [21], see also Refs.[18][19][20] for the cases of certain small n.However, kernel functions K n where previously only known approximately, for situations involving only a low order of anomalous terms, see the discussion in Sec. 5. We define an anomalous term of order a = 1, 2, ...n − 1 as a summand contributing to K n (Ω 1 , ..., Ω n−1 ) that contains a product of a Kronecker delta functions δ 0,x , where x is a sum of a subset of {Ω 1 , ..., Ω n−1 }.As can be seen in Fig. 1(c), these anomalous contributions to G A 1 ...A n ω 1 , ..., ω n−1 correspond to qualitatively important sharp features.
In the next section, we present a simple, exact expression for general K n (Ω 1 , ..., Ω n−1 ).Readers not interested in the derivation can directly skip to the result in Eq. ( 26) or its explicit form for n = 2, 3, 4, 5 in Sec. 5.

General kernel function
Assuming the spectrum and matrix elements entering Eq. ( 15) are known, the remaining task is to find expressions for the kernel function K n (Ω 1 , ..., Ω n−1 ) defined via Eqns.( 13) and ( 14) as the part of K n (Ω 1 , Ω 2 , ..., Ω n ) multiplying βδ 0,Ω 1 +Ω 2 +...+Ω n .To facilitate the presentation in this section, in Eq. ( 13) we rename the integration variables τ k → τ n−k+1 and define new arguments z n− j+1 = Ω j for j = 1, 2, ..., n − 1, As indicated in Eq. ( 16), we call h k (τ k ) the integrand for the integrand is given by h 1 (τ 1 ) = e z 1 τ 1 and we will find h k for k = 2, 3, ..., n iteratively.For z ∈ , we define the abbreviations δ z ≡ δ 0,z and and consider the integral (for p = 0, 1, 2, ... and τ ≥ 0, proof by partial integration and induction) Recall that we are only interested in the contribution K n (z n , z n−1 , ..., z 2 ) that fulfills frequency conservation, see Eq. ( 17).The δ z 1 +z 2 +...+z n in front of this term arises from the final τ n integration of h n (τ n ) ∝ e (z 1 +z 2 +...+z n )τ n via the first term in Eq. ( 19).This however requires that all z k (except the vanishing ones, of course) remain in the exponent during the iterative integrations.This requirement is violated by the last term in the general integral (19) (which comes from the lower boundary of the integral).All terms in K n that stem from this last term in Eq. ( 19) thus contribute to R n and can be dropped in the following [21].Note however, that it is straightforward to generalize our approach and keep these terms if the full K n is required.
To define the iterative procedure to solve the n-fold integral in Eq. ( 16), we make the ansatz which follows from the form of the integral (19) and our decision to disregard the terms contributing to R n .The ansatz (20) is parameterized by the numbers f k (l) with l = 0, 1, ..., k − 1.
These numbers have to be determined iteratively, starting from f k=1 (l = 0) = 1, read off from h 1 (τ 1 ) = e z 1 τ 1 , c.f. Eq. ( 16).Iteration rules to obtain the f k (l) from f k−1 (l) are easily derived from Eqns. ( 16), ( 19) and (20).We obtain the recursion relation This can be understood as a matrix-vector product of where ∆k ≡ ∆ z k +...+z 2 +z 1 , δk ≡ δ z k +...+z 2 +z 1 .The tilde on top of the δk and ∆k signals the presence of a sum of z j in the arguments (below we will define related quantities without tilde for the sum of Ω j ).Note that the first (second) term in brackets of Eq. ( 22) comes from the first (second) term in square brackets of Eq. (19).
The next step is to find K n (z n , z n−1 , ..., z 2 ).This requires to do the integral β 0 dτ n h n (τ n ) which can be again expressed via Eq.( 19) but with the replacement τ → β.Only the first term provides a βδ z 1 +z 2 +...+z n and is thus identified with K n .We find: The argument z 1 that the right hand side of Eq. ( 23) depends on is to be replaced by . Then, to conform with Eq. ( 15), we reinstate Ω j = z n− j+1 for j = 1, 2, ..., n − 1.This amounts to replacing the terms δj and ∆j that appear in f n (l) as follows, where we used Finally, we can express Eq. ( 23) using a product of n − 1 matrices M multiplying the initial length-1 vector with entry f 1 (0) = 1.
Transferring to the Ω-notation by using Eqns.( 24) and (25), we obtain The closed form expression (26) of the universal kernel, to be used in the spectral representation (15), is our main result.By definition it is free of any singularities as the case of vanishing denominators is explicitly excluded in Eq. ( 18).
5 Explicit kernel functions K n Ω 1 , ..., Ω n−1 for n = 2, 3, 4, 5 While the previous section gives a closed form expression for kernel functions of arbitrary order, we here evaluate the universal kernel functions K n (Ω 1 , ..., Ω n−1 ) defined in Eq. ( 14) from Eq. ( 26) for n = 2, 3, 4, 5 and show the results in Tab. 1.In each column, the kernel function in the top row is obtained by first multiplying the entries listed below it in the same column by the common factor in the rightmost column and then taking the sum.The symbols δ j and ∆ j for j = 1, 2, ..., n − 1 which appear in Tab. 1 are defined by compare also to the previous section.As an example, for n = 2 and n = 3 we obtain from Tab. 1 respectively.The rows of Tab. 1 are organized with respect to the number a of factors δ l in the summands.Here, a = 0 indicates the regular part and a = 1, 2, ..., n − 1 indicates anomalous terms.There are n − 1 choose a anomalous terms of order a.Our results are exact and go substantially beyond existing expressions in the literature -these are limited to n ≤ 3 [18] or to fermionic n = 4 [19][20][21] with a = 0, 1 (and a = 2, 3 guaranteed to vanish, see below) or arbitrary n with a = 0 [21].Alternative expressions for the n = 3, 4 kernel functions with a ≤ 1 were given in [21], but they are consistent with our kernel functions as they yield the same correlation functions, see the Appendix.
In the case of purely fermionic correlators (all A k fermionic), individual Matsubara frequencies ω k cannot be zero and thus the Ω 12) always have a finite imaginary part and are non-zero, regardless of the eigenenergies.In this case, only sums of an even number of frequencies can be zero, and we can simplify δ 1 = δ 3 = δ 5 = ... = 0.The expressions for the kernels in Tab. 1, now denoted by K n | F for the fermionic case, simplify to Table 1: Universal kernel functions K n Ω 1 , ..., Ω n−1 for n = 2, 3, 4, 5 defined in Eq. ( 14) and calculated from Eq. ( 26) in Sec. 4. In each column, the kernel function in the top row is obtained by first multiplying the entries listed below it in the same column by the common factor in the rightmost column and then taking the sum, see Eqns. ( 30) and ( 31) as examples.The symbols δ j and ∆ j appearing are defined in Eqns.( 28) and ( 29).The rows are organized with respect to the number a of appearances of δ j , i.e. the order of the anomalous terms.
This concludes the general part of this work.Next, we consider two example systems frequently discussed in the condensed matter theory literature.Using our formalism, we provide analytical forms of correlation functions that to the best of our knowledge were not available before.
6 Applications: Hubbard atom and free spin

Fermionic Hubbard atom
The Hubbard atom (HA) describes an isolated impurity or otherwise localized system with Hamiltonian see Fig. 1(b) for a sketch.The HA corresponds to the limit of vanishing system-bath coupling of the Anderson impurity model (AIM), or vanishing hopping in the Hubbard model (HM).The particle number operators n σ = d † σ d σ count the number of fermionic particles with spin σ ∈ {↑, ↓}, each contributing an onsite energy ε shifted by an external magnetic field h in z-direction.An interaction energy U is associated to double occupation.
Due to its simplicity and the four-dimensional Hilbert space, the correlation functions for the HA can be found analytically using the spectral representation.It is therefore often used for benchmarking [3,27,28].The presence of the interaction term leads to a non-vanishing n = 4 one-line irreducible vertex function.The HA serves as an important reference point to study and interpret properties of the AIM and HM beyond the one-particle level, for example diver-gences of two-line irreducible vertex functions [29][30][31][32] and signatures of the local moment formation in generalized susceptibilities [16,33].Using the fermionic kernels in Eqns.( 32) and ( 33), we have checked that our formalism reproduces the results for the 2-point and 4point correlators given in Refs.[19,21,26] for half-filling, ε = −U/2 and h = 0.
Correlation functions including bosonic operators describe the asymptotic behaviour of the n = 4 fermion vertex for large frequencies [34] or the interaction of electrons by the exchange of an effective boson [35,36].These relations involve correlation functions of two bosonic operators or of one bosonic and two fermionic operators, giving rise to expressions possibly anomalous in at most one frequency argument, i.e. a ≤ 1.
For the HA, AIM and HM, bosonic correlation functions for n > 2 have not been considered thoroughly so far.Only recently, steps in this direction were taken, particularly in the context of non-linear response theory [3].The response of a system to first and second order in an external perturbation is described by 2-and 3-point correlation functions, respectively.For the HA, physically motivated perturbations affect the onsite energy via a term δ ε n or take the form of a magnetic field δ h • S. Here, the parameters δ ε and δ h denote the strength of the perturbation and we define The resulting changes of the expectation values of the density or magnetization in arbitrary direction are described in second order of the perturbation by connected parts of the correlation functions G A 1 A 2 A 3 (τ 1 , τ 2 , τ 3 ), with A i ∈ {n, S x , S y , S z }, where the time-ordered expectation value is evaluated with respect to the unperturbed system (35) and Fourier transformed to the frequencies of interest.These objects have been studied numerically in Ref. [3].In the following, we give explicit, analytic expressions of the full correlation functions G A 1 A 2 A 3 (ω 1 , ω 2 ) (i.e.including disconnected parts), for arbitrary parameters ε, U and h and for all possible operator combinations using the (bosonic) kernel function K 3 , see Eq. (31).To the best of our knowledge, these expressions have not been reported before.
The eigenstates of the HA Hamiltonian (35) [see Fig.
and obtain all non-vanishing bosonic 3-point correlation functions (where We observe that each conserved quantity, in this case n and S z , contributes an anomalous term ∝ δ ω k in its respective frequency argument ω k .If an operator A k is conserved [H, A k ] = 0, the basis over which we sum in Eq. ( 15) can be chosen such that both H and A k are diagonal, A k ̸ = 0 for some state |1〉 the vanishing eigenenergy difference leads to the appearance of an anomalous contribution.If the operators in the correlator additionally commute with each other, in our case for example [n, S z ] = 0, there exists a basis in which all operators and the Hamiltonian are diagonal, giving rise to correlation functions anomalous in all frequency arguments.
In the limit of vanishing field h → 0, we introduce an additional degeneracy E ↑ = E ↓ = ε in the system, potentially resulting in additional anomalous contributions.The corresponding correlation functions can then be obtained in two ways.Either we recompute them using the kernel function K 3 or we take appropriate limits, for example resulting in with all other correlation functions vanishing.As already pointed out in Ref. [3], only the last correlation function retains a nontrivial frequency dependence due to non-commuting operators.

Free spin S
We now consider correlation functions of a free spin of length S, without a magnetic field, so that temperature T = 1/β is the only energy scale.The operators {S α } α=x, y,z fulfill S x S x + S y S y + S z S z = S(S + 1) and the SU(2) algebra [S α 1 , S α 2 ] = i α 3 ={x, y,z} ε α 1 α 2 α 3 S α 3 , thus they are bosonic.Since the Hamiltonian vanishes and therefore all eigenenergies are zero, every Ω a b k in the spectral representation (15) can vanish and a proper treatment of all anomalous terms is essential.As the Heisenberg time dependence is trivial, S α (τ) = S α , the non-trivial frequency dependence of the correlators, which can be non-vanishing at any order n > 1, derives solely from the action of imaginary time-ordering.
The correlators are required, for example, as the non-trivial initial condition for the spin-fRG recently suggested by Kopietz et al.,Refs. [13,[37][38][39][40].However, for n > 3 they are so far only partially available: They are either given for restricted frequency combinations, or for the purely classical case S α 1 = S α 2 = ... = S α n where the SU(2) algebra does not matter, or for finite magnetic field via an equation of motion [37] or diagrammatic approach [41,42].literature. 1 We also confirmed the classical result for G S z S z S z S z , which in our full quantum formalism requires some non-trivial cancellations.To arrive at our results, we used the identity We finally comment on the relation between the n = 3 free spin-S correlator G S + S − S z from Tab.
(2) and the result for G S x S y S z found for the zero-field limit of the HA in Eq. ( 49).The operators S x, y,z for the Hubbard model [c.f.Eq. ( 36)] project to the singly-occupied S = 1/2 subspace spanned by the states ↑ , ↓ .Thus, using G S x S y S z = iG S + S − S z and specializing the free spin result from Tab. (2) to S = 1/2 (where b 1 = 1/4) we find agreement with the HA result (49) up to the factor 2e −βε /Z.This factor represents the expectation value of the projector to the singly-occupied sector in the HA Hilbert space and goes to unity in the localmoment regime.

Conclusion
In summary, we have provided exact universal kernel functions for the spectral representation of the n-point Matsubara correlator.Our results are an efficient alternative to equation-ofmotion approaches which often have difficulties to capture anomalous terms related to conserved or commuting operators.We expect our results to be useful for various benchmarking applications, as starting points for emerging many-body methods and for unraveling the physical interpretation of n-point functions in various settings.Our results also apply in the limit T → 0 where the formally divergent anomalous contributions are to be understood as βδ ω,0 → 2πδ(ω).Some of these Dirac delta-functions will vanish after subtracting the disconnected contributions, others indicate truely divergent susceptibilities like the 1/T Curie law for the spin-susceptiblity of the Hubbard atom in the local moment regime [26].Although our work has focused on imaginary frequency (Matsubara) correlators, with analytical expressions now at hand, it is also interesting to study the intricacies of analytical continuation to real frequencies and thus to further explore the connection of Matsubara and Keldysh correlators [43].
k H denotes Heisenberg time evolution.Here and in the following, k = 1, 2, ..., n.The expectation value is calculated as ... = tr[ρ...] where ρ = exp(−β H)/Z is the thermal density operator at inverse temperature β = 1/T and Z = tr exp(−β H) is the partition function.Note that other conventions for the n-point function differing by a prefactor are also used in the literature, e.g.Ref.

7 )
the the second term of p = 123 and the first term of p = 231 as well as the second term of p = 231 and the first term of p = 312 cancel, over the second set of cyclically related permutations p = 132, 213, 321 leads to a vanishing result, leading to the conclusion that 1 Thus we have shown that both kernel functions in Eqns.(A.1) and (A.2) are equivalent as they yield the same correlation functions after summing over all permutations.The same statement holds true for case of n = 4 and a = 1.