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

Johannes Halbinger, Benedikt Schneider and Björn Sbierski  Department of
Physics and Arnold Sommerfeld Center for Theoretical Physics (ASC), Ludwig-
Maximilians-Universität München, Theresienstr. 37, München D-80333, Germany
Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4,
D-80799 München, Germany

(August 28, 2024)

###### Abstract

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.

##  I Introduction and definitions

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 (Negele and Orland, 1988). 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
angle-resolved photon emission spectroscopy (Bruus and Flensberg, 2004). 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 (Kappl et al.,
2023). In the emerging field of cold atomic quantum simulation, (equal-time)
$n$-point functions are even directly accessible (Semeghini et al., 2021).

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$ (Hedin, 1965; Bickers, ).
Although these relations can usually not be solved exactly, they form a
valuable starting point for further methodological developments like the
parquet approximation (Roulet et al., 1969). Thus even if the 4-point
correlator (or, in that context, its essential part, the one-line irreducible
vertex (Negele and Orland, 1988)) 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 (Kopietz et al., 2010; Metzner et al., 2012). 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 (Wetterich, 1993). 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 (Negele and Orland, 1988). 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 $n$-point 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 (Lee et al., 2021).

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 (Reuther and
Thomale, 2014; Krieg and Kopietz, 2019; Rückriegel et al., 2022), or appear
intrinsically in a method like in diagrammatic extensions of dynamical mean
field theory (Rohringer et al., 2018) 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. (Chalupa et al., 2021)
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 (Lehmann, 1954; Negele
and Orland, 1988) 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. (Shvaika, 2006; Hafermann et al., 2009; Shvaika,
2016) and in particular the recent approach by Kugler _et al._ (Kugler et al.,
2021), 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. (Kugler et
al., 2021)). 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 (Shvaika,
2006), for the fermionic 4-point case (Hafermann et al., 2009; Shvaika, 2016;
Kugler et al., 2021) or for the general $n$-point case (Kugler et al., 2021)
but disregarding anomalous contributions to the sum that the kernel function
consists 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 long-term 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 (Kwok and
Schultz, 1969; Watzenböck et al., 2022).

The structure of the paper is as follows: In Sec. II we define the Matsubara
$n$-point function
$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n-1}\bigr{)}$ and review some
of its properties. The spectral representation is derived in Sec. III with Eq.
(15) being the central equation written in terms of the kernel function
$K_{n}(\Omega_{1},...,\Omega_{n-1})$. Our main result is an exact closed-form
expression of this most general kernel function which is given in Sec. IV.
Examples for $n=2,3,4,5$ are given in Sec. V 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. VI, 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. VII.

Figure 1: (a) Ordering convention for imaginary times in Eq. (9). (b)
Eigenstates and energies of the Hubbard atom. (c) Matsubara correlation
function $G_{S^{x}S^{y}S^{z}}(\omega_{1},\omega_{2})$ with $\omega_{j}=2\pi
m_{j}/\beta$ ($m_{j}\in\mathbb{Z}$, $j=1,2$) for the Hubbard atom (35) at
$\beta=10$, $h=0.1$, $\epsilon=-2$, $U=2$, see Eq. (45). The sharp anti-
diagonal ray $\propto\delta_{\omega_{1}+\omega_{2},0}$ represents an anomalous
term of order $a=1$. The other broadened rays become sharp and anomalous for
$h\rightarrow 0$, see Eq. (49).

##  II Definition of Matsubara $n$-point function
$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n}\bigr{)}$

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, $A_{1}=d^{\dagger}d\equiv n$, $A_{2}=d$, $A_{3}=d^{\dagger}$ where
$d^{\dagger}$ 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. (Tsvelik, 2007).
Spin operators are thus bosonic.

We define the imaginary time-ordered $n$-point correlation functions for
$\tau_{k}\in[0,\beta]$ (Rohringer et al., 2012; Rohringer, 2013),

$G_{A_{1}A_{2}...A_{n}}\left(\tau_{1},\tau_{2},...,\tau_{n}\right)\equiv\left\langle\mathcal{T}A_{1}(\tau_{1})A_{2}(\tau_{2})...A_{n}(\tau_{n})\right\rangle,$
(1)

where $A_{k}(\tau_{k})=e^{\tau_{k}H}A_{k}e^{-\tau_{k}H}$ denotes Heisenberg
time evolution. Here and in the following, $k=1,2,...,n$. The expectation
value is calculated as $\bigl{\langle}...\bigr{\rangle}=\mathrm{tr}[\rho...]$
where $\rho=\exp(-\beta H)/Z$ is the thermal density operator at temperature
$\beta=1/T$ and $Z=\mathrm{tr}\exp(-\beta H)$ denotes 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. (Kugler et al., 2021)
multiplies with $(-1)^{n-1}$. In Eq. (1), the imaginary time-ordering operator
$\mathcal{T}$ orders the string of Heisenberg operators,

$\mathcal{T}A_{1}(\tau_{1})A_{2}(\tau_{2})...A_{n}(\tau_{n})\equiv\boldsymbol{\zeta}(p)A_{p(1)}(\tau_{p(1)})A_{p(2)}(\tau_{p(2)})...A_{p(n)}(\tau_{p(n)}),$
(2)

where $p$ is the permutation $p\in S_{n}$ such that
$\tau_{p(1)}>\tau_{p(2)}>...>\tau_{p(n)}$ [see Fig. 1(a)] and the sign
$\boldsymbol{\zeta}(p)$ is $-1$ if the operator string
$A_{p(1)}A_{p(2)}...A_{p(n)}$ differs from $A_{1}A_{2}...A_{n}$ by an odd
number of transpositions of fermionic operators, otherwise it is $+1$. The
special case $n=2$, with $\boldsymbol{\zeta}(12)=1$ and
$\boldsymbol{\zeta}(21)=\zeta$ ($\zeta=1$ for $A_{1,2}$ bosonic, $\zeta=-1$
for $A_{1,2}$ fermionic), simplifies to

$\mathcal{T}A_{1}(\tau_{1})A_{2}(\tau_{2})=\begin{cases}A_{1}(\tau_{1})A_{2}(\tau_{2})&:\tau_{1}>\tau_{2},\\\
\zeta A_{2}(\tau_{2})A_{1}(\tau_{1})&:\tau_{2}>\tau_{1}.\end{cases}$ (3)

Imaginary time-ordered correlation functions (1) fulfill certain properties
which we review in the following, see e.g. (Rohringer, 2013) for a more
extensive discussion. First, they are invariant under translation of all time
arguments,

$G_{A_{1}A_{2}...A_{n}}\left(\tau_{1},\tau_{2},...,\tau_{n}\right)=G_{A_{1}A_{2}...A_{n}}\left(\tau_{1}+\tau,\tau_{2}+\tau,...,\tau_{n}+\tau\right),$
(4)

with $\tau\in\mathbb{R}$ such that $\tau_{k}+\tau\in[0,\beta]$. They also
fulfill periodic or anti-periodic boundary conditions for the individual
arguments $\tau_{k}$,

$G_{A_{1}...A_{n}}\left(\tau_{1},...,\tau_{k}=0,...,\tau_{n}\right)=\zeta_{k}G_{A_{1}...A_{n}}\left(\tau_{1},...,\tau_{k}=\beta,...,\tau_{n}\right)$
(5)

where $\zeta_{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,

$\displaystyle G_{A_{1}...A_{n}}\left(\tau_{1},...,\tau_{n}\right)$
$\displaystyle\equiv$
$\displaystyle\beta^{-n}\sum_{\omega_{1},...,\omega_{n}}e^{-i(\omega_{1}\tau_{1}+...+\omega_{n}\tau_{n})}G_{A_{1}...A_{n}}\left(\omega_{1},...,\omega_{n}\right),$
(6) $\displaystyle G_{A_{1}...A_{n}}\left(\omega_{1},...,\omega_{n}\right)$
$\displaystyle=$
$\displaystyle\int_{0}^{\beta}\mathrm{d}\tau_{1}\cdots\int_{0}^{\beta}\mathrm{d}\tau_{n}e^{+i(\omega_{1}\tau_{1}+...+\omega_{n}\tau_{n})}G_{A_{1}...A_{n}}\left(\tau_{1},...,\tau_{n}\right),$
(7)

where $\omega_{k}=2\pi m_{k}/\beta$ or $\omega_{k}=2\pi(m_{k}+1/2)/\beta$ with
$m_{k}\in\mathbb{Z}$ bosonic or fermionic Matsubara frequencies, respectively,
and $\sum_{\omega_{k}}$ is shorthand for $\sum_{m_{k}\in\mathbb{Z}}$. 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}}(\omega_{1},...,\omega_{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),

$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n-1},\omega_{n}\bigr{)}\equiv\beta\delta_{0,\omega_{1}+...+\omega_{n}}G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n-1}\bigr{)},$
(8)

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 $\beta$ and the Kronecker delta function.

##  III Spectral representation of
$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n-1}\bigr{)}$

The integrals involved in the Fourier transformation (7) generate all $n!$
different orderings of the time arguments $\tau_{k}$. As in Ref. (Kugler et
al., 2021) it is thus convenient to use a sum over all $n!$ permutations $p\in
S_{n}$ and employ a product of $n-1$ step-functions $\theta$, with
$\theta(x)=1$ for $x>0$ and $0$ otherwise, to filter out the unique ordering
for which $\beta>\tau_{p(1)}>\tau_{p(2)}>...>\tau_{p(n-1)}>\tau_{p(n)}>0$, see
Fig. 1(a),

$\displaystyle G_{A_{1}...A_{n}}(\tau_{1},...,\tau_{n})$ $\displaystyle=$
$\displaystyle\sum_{p\in
S_{n}}\boldsymbol{\zeta}(p)\left[\prod_{i=1}^{n-1}\theta(\tau_{p(i)}-\tau_{p(i+1)})\right]\left\langle
A_{p(1)}(\tau_{p(1)})A_{p(2)}(\tau_{p(2)})...A_{p(n)}(\tau_{p(n)})\right\rangle.$
(9)

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
$\bigl{|}j_{1}\bigr{\rangle},\bigl{|}j_{2}\bigr{\rangle},...$ and
$E_{j_{1}},E_{j_{2}},...$ we employ
$\bigl{|}\underline{1}\bigr{\rangle},\bigl{|}\underline{2}\bigr{\rangle},...$
and $E_{\underline{1}},E_{\underline{2}}$,…. for compressed notation and
denote operator matrix elements as
$A^{\underline{1}\underline{2}}=\left\langle\underline{1}|A|\underline{2}\right\rangle$.
We obtain

$\displaystyle G_{A_{1}...A_{n}}(\tau_{1},...,\tau_{n})=\sum_{p\in
S_{n}}\boldsymbol{\zeta}(p)\left[\prod_{i=1}^{n-1}\theta(\tau_{p(i)}-\tau_{p(i+1)})\right]$
(10)
$\displaystyle\times\frac{1}{Z}\sum_{\underline{1}...\underline{n}}e^{-\beta
E_{\underline{1}}}e^{\tau_{p(1)}E_{\underline{1}}}A_{p(1)}^{\underline{1}\,\underline{2}}e^{(-\tau_{p(1)}+\tau_{p(2)})E_{\underline{2}}}A_{p(2)}^{\underline{2}\,\underline{3}}e^{(-\tau_{p(2)}+\tau_{p(3)})E_{\underline{3}}}...e^{(-\tau_{p(n-1)}+\tau_{p(n)})E_{\underline{n}}}A_{p(n)}^{\underline{n}\,\underline{1}}e^{-\tau_{p(n)}E_{\underline{1}}},$

and apply the Fourier transform according to the definition (7),

$\displaystyle
G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n}\bigr{)}=\frac{1}{Z}\sum_{p\in
S_{n}}\boldsymbol{\zeta}(p)\sum_{\underline{1}...\underline{n}}e^{-\beta
E_{\underline{1}}}A_{p(1)}^{\underline{1}\,\underline{2}}A_{p(2)}^{\underline{2}\,\underline{3}}...A_{p(n)}^{\underline{n}\,\underline{1}}$
(11)
$\displaystyle\times\left[\int_{0}^{\beta}\\!\\!\mathrm{d}\tau_{p(1)}e^{\Omega_{p(1)}^{\underline{1}\,\underline{2}}\tau_{p(1)}}\right]\left[\int_{0}^{\tau_{p(1)}}\\!\\!\mathrm{d}\tau_{p(2)}e^{\Omega_{p(2)}^{\underline{2}\,\underline{3}}\tau_{p(2)}}\right]...\left[\int_{0}^{\tau_{p(n-2)}}\\!\\!\mathrm{d}\tau_{p(n-1)}e^{\Omega_{p(n-1)}^{\underline{n-1}\,\underline{n}}\tau_{p(n-1)}}\right]\left[\int_{0}^{\tau_{p(n-1)}}\\!\\!\mathrm{d}\tau_{p(n)}e^{\Omega_{p(n)}^{\underline{n}\,\underline{1}}\tau_{p(n)}}\right],$

where we defined

$\Omega_{k}^{\underline{a}\,\underline{b}}\equiv
i\omega_{k}+E_{\underline{a}}-E_{\underline{b}}\in\mathbb{C}.$ (12)

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 second line can be
regarded as a universal kernel function defined for general
$\\{\Omega_{1},\Omega_{2},...,\Omega_{n}\\}$ probed at
$\Omega_{k}\in\mathbb{C}$ which depends on the system and correlators via
(12). Upon renaming the $\tau$-integration variables
$\tau_{p(k)}\rightarrow\tau_{k}$, this kernel function is written as follows:

$\displaystyle\mathcal{K}_{n}\left(\Omega_{1},...,\Omega_{n}\right)$
$\displaystyle\equiv$
$\displaystyle\left[\int_{0}^{\beta}\mathrm{d}\tau_{1}e^{\Omega_{1}\tau_{1}}\right]\left[\int_{0}^{\tau_{1}}\mathrm{d}\tau_{2}e^{\Omega_{2}\tau_{2}}\right]...\left[\int_{0}^{\tau_{n-2}}\mathrm{d}\tau_{n-1}e^{\Omega_{n-1}\tau_{n-1}}\right]\left[\int_{0}^{\tau_{n-1}}\mathrm{d}\tau_{n}e^{\Omega_{n}\tau_{n}}\right]$
(13) $\displaystyle\equiv$
$\displaystyle\beta\delta_{0,\Omega_{1}+\Omega_{2}+...+\Omega_{n}}K_{n}\left(\Omega_{1},...,\Omega_{n-1}\right)+R_{n}\left(\Omega_{1},...,\Omega_{n}\right).$
(14)

In the second line we split $\mathcal{K}_{n}$ into a part $K_{n}$ proportional
to $\beta\delta_{0,\Omega_{1}+\Omega_{2}+...+\Omega_{n}}$ and the rest
$R_{n}$. We dropped $\Omega_{n}$ from the argument list of $K_{n}$ which can
be reconstructed from $\\{\Omega_{1},...,\Omega_{n-1}\\}$.

Finally, we express
$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n}\bigr{)}$ of Eq. (11) using
the kernel $\mathcal{K}_{n}$ so that the general $\Omega_{k}\in\mathbb{C}$ get
replaced by $\Omega_{k}^{\underline{a}\,\underline{b}}$ of Eq. (12). For
these,
$\Omega_{p(1)}^{\underline{1}\,\underline{2}}+\Omega_{p(2)}^{\underline{2}\,\underline{3}}+...+\Omega_{p(n)}^{\underline{n}\,\underline{1}}=i(\omega_{1}+\omega_{2}+...+\omega_{n})$,
since the $E_{\underline{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\in S_{n}$, so that only the terms proportional to $K_{n}$ remain. We drop
the $\beta\delta_{0,\omega_{1}+\omega_{2}+...+\omega_{n}}$ from both sides
[c.f. Eq. (8)] and find the spectral representation of the $n$-point
correlation function in the Matsubara formalism,

$\boxed{G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n-1}\bigr{)}=\frac{1}{Z}\sum_{p\in
S_{n}}\boldsymbol{\zeta}(p)\sum_{\underline{1}...\underline{n}}e^{-\beta
E_{\underline{1}}}A_{p(1)}^{\underline{1}\,\underline{2}}A_{p(2)}^{\underline{2}\,\underline{3}}...A_{p(n)}^{\underline{n}\,\underline{1}}\times
K_{n}\left(\Omega_{p(1)}^{\underline{1}\,\underline{2}},\Omega_{p(2)}^{\underline{2}\,\underline{3}},...,\Omega_{p(n-1)}^{\underline{n-1}\,\underline{n}}\right).}$
(15)

An equivalent expression was derived in the literature before (Kugler et al.,
2021), see also Refs. (Shvaika, 2006; Hafermann et al., 2009; Shvaika, 2016)
for the cases of certain small $n$. However, the kernel functions $K_{n}$
where previously only known approximately, for situations involving only a low
order of anomalous terms, see the discussion in Sec. V. We define an
_anomalous_ term of order $a=1,2,...n-1$ as a summand contributing to
$K_{n}\left(\Omega_{1},...,\Omega_{n-1}\right)$ that contains a product of $a$
Kronecker delta functions $\delta_{0,x}$, where $x$ is a sum of a subset of
$\\{\Omega_{1},...,\Omega_{n-1}\\}$. As can be seen in Fig. 1(c), these
anomalous contributions to
$G_{A_{1}...A_{n}}\bigl{(}\omega_{1},...,\omega_{n}\bigr{)}$ correspond to
qualitatively important sharp features.

In the next section, we present a simple, exact expression for general
$K_{n}\left(\Omega_{1},...,\Omega_{n-1}\right)$. 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. (V).

##  IV General kernel function
$K_{n}\bigl{(}\Omega_{1},...,\Omega_{n-1}\bigr{)}$

Assuming the spectrum and matrix elements entering Eq. (15) are known, the
remaining task is to find expressions for the kernel function
$K_{n}\left(\Omega_{1},...,\Omega_{n-1}\right)$ defined via Eqns. (13) and
(14) as the part of
$\mathcal{K}_{n}\left(\Omega_{1},\Omega_{2},...,\Omega_{n}\right)$ multiplying
$\beta\delta_{0,\Omega_{1}+\Omega_{2}+...+\Omega_{n}}$. To facilitate the
presentation in this section, in Eq. (13) we rename the integration variables
$\tau_{k}\rightarrow\tau_{n-k+1}$ and define new arguments
$z_{n-j+1}=\Omega_{j}$ for $j=1,2,...,n-1$,

$\displaystyle\mathcal{K}_{n}\left(\Omega_{1}=z_{n},\Omega_{2}=z_{n-1},...,\Omega_{n}=z_{1}\right)$
$\displaystyle=\Bigl{[}\int_{0}^{\beta}\mathrm{d}\tau_{n}e^{z_{n}\tau_{n}}\Bigr{]}\Bigl{[}\int_{0}^{\tau_{n}}\mathrm{d}\tau_{n-1}e^{z_{n-1}\tau_{n-1}}\Bigr{]}...\Bigl{[}\int_{0}^{\tau_{3}}\mathrm{d}\tau_{2}\underset{\equiv
h_{2}(\tau_{2})}{\underbrace{e^{z_{2}\tau_{2}}\Bigr{]}\Bigl{[}\int_{0}^{\tau_{2}}\mathrm{d}\tau_{1}\underset{\equiv
h_{1}(\tau_{1})}{\underbrace{e^{z_{1}\tau_{1}}}}\Bigr{]}}}$ (16)
$\displaystyle=\beta\delta_{0,z_{1}+z_{2}+...+z_{n}}K_{n}\left(z_{n},z_{n-1},...,z_{2}\right)+R\left(z_{n},z_{n-1},...,z_{1}\right).$
(17)

As indicated in Eq. (16), we call $h_{k}(\tau_{k})$ the integrand of the
$\int_{0}^{\tau_{k+1}}\mathrm{d}\tau_{k}$ integral for $k=1,2,...,n$. At $k=1$
this integrand is given by $h_{1}(\tau_{1})=e^{z_{1}\tau_{1}}$ and we will
find $h_{k}$ for $k=2,3,...,n$ iteratively. For $z\in\mathbb{C}$, we define
the abbreviations $\delta_{z}\equiv\delta_{0,z}$ and

$\Delta_{z}\equiv\begin{cases}0&\textnormal{{if} }z=0\\\
\frac{{1}}{z}&\textnormal{{if} }z\neq 0\end{cases}$ (18)

and consider the integral (for $p=0,1,2,...$ and $\tilde{\tau}\geq 0$, proof
by partial integration and induction)

$\int_{0}^{\tilde{\tau}}\mathrm{d}\tau\,\tau^{p}e^{z\tau}=\left[\frac{\tilde{\tau}^{p+1}}{p+1}\delta_{z}+p!\left(-1\right)^{p}\Delta_{z}^{1+p}\sum_{l=0}^{p}\frac{(-1)^{l}}{l!}\Delta_{z}^{-l}\tilde{\tau}^{l}\right]e^{z\tilde{\tau}}-p!\left(-1\right)^{p}\Delta_{z}^{p+1}.$
(19)

Recall that we are only interested in the contribution
$K_{n}\left(z_{n},z_{n-1},...,z_{2}\right)$ that fulfills frequency
conservation, see Eq. (17). The $\delta_{z_{1}+z_{2}+...+z_{n}}$ in front of
this term arises from the final $\tau_{n}$ integration of
$h_{n}(\tau_{n})\propto e^{(z_{1}+z_{2}+...+z_{n})\tau_{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
$\mathcal{K}_{n}$ that stem from this last term in Eq. (19) thus contribute to
$R_{n}$ and can be dropped in the following (Kugler et al., 2021). Note
however, that it is straightforward to generalize our approach and keep these
terms if the full $\mathcal{K}_{n}$ is required.

To define the iterative procedure to solve the $n$-fold integral in Eq. (16),
we make the ansatz

$h_{k}(\tau_{k})=\sum_{l=0}^{k-1}f_{k}(l)\tau_{k}^{l}e^{\left(z_{k}+z_{k-1}+...+z_{1}\right)\tau_{k}},$
(20)

which follows from the form of the integral (19) and our decision to disregard
the terms contributing to $R_{n}$. The ansatz $\eqref{eq:h_k-Ansatz}$ 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}(\tau_{1})=e^{z_{1}\tau_{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 $\eqref{eq:h_k-Ansatz}$. We obtain the recursion relation

$f_{k}(l)=\sum_{p=0}^{k-1}\tilde{M}_{k-1}(l,p)f_{k-1}(p)$ (21)

written as a matrix-vector product of
$\mathbf{f}_{k-1}=(f_{k-1}(0),f_{k-1}(1),...,f_{k-1}(k-2))^{\mathrm{T}}$ with
the $k\times(k-1)$-matrix

$\tilde{M}_{k-1}(l,p)=\frac{p!}{l!}\left[\delta_{l,p+1}\tilde{\delta}_{k-1}+\theta\left(p-l+1/2\right)\left(-1\right)^{l+p}\tilde{\Delta}_{k-1}^{1+p-l}\right],$
(22)

where $\tilde{\Delta}_{k}\equiv\Delta_{z_{k}+...+z_{2}+z_{1}}$,
$\tilde{\delta}_{k}\equiv\delta_{z_{k}+...+z_{2}+z_{1}}$. The tilde on top of
the $\tilde{\delta}_{k}$ and $\tilde{\Delta}_{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 $\Omega_{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}\left(z_{n},z_{n-1},...,z_{2}\right)$. This
requires to do the integral
$\int_{0}^{\beta}\mathrm{d}\tau_{n}h_{n}(\tau_{n})$ which can be again
expressed via Eq. (19) but with the replacement
$\tilde{\tau}\rightarrow\beta$. Only the first term provides a
$\beta\delta_{z_{1}+z_{2}+...+z_{n}}$ and is thus identified with $K_{n}$. We
find:

$K_{n}\left(z_{n},z_{n-1},...,z_{2}\right)=\sum_{l=0}^{n-1}\frac{\beta^{l}f_{n}(l)}{l+1}.$
(23)

The argument $z_{1}$ that the right hand side of Eq. (23) depends on is to be
replaced by $z_{1}=-z_{2}-z_{3}-...-z_{n}$, in line with the arguments in
$K_{n}\left(z_{n},z_{n-1},...,z_{2}\right)$. Then, to conform with Eq. (15),
we reinstate $\Omega_{j}=z_{n-j+1}$ for $j=1,2,...,n-1$. This amounts to
replacing the terms $\tilde{\delta}_{j}$ and $\tilde{\Delta}_{j}$ that appear
in $f_{n}(l)$ as follows,

$\displaystyle\tilde{\delta}_{j}$ $\displaystyle=$
$\displaystyle\delta_{z_{j}+...+z_{2}+z_{1}}=\delta_{\Omega_{1}+\Omega_{2}+...+\Omega_{n-j}}\equiv\delta_{n-j},$
(24) $\displaystyle-\tilde{\Delta}_{j}$ $\displaystyle=$
$\displaystyle-\Delta_{z_{j}+...+z_{2}+z_{1}}=\Delta_{\Omega_{1}+\Omega_{2}+...+\Omega_{n-j}}\equiv\Delta_{n-j},$
(25)

where we used $\Omega_{1}+\Omega_{2}+...+\Omega_{n}=0=z_{n}+...+z_{2}+z_{1}$.
Finally, we can express Eq. (23) using a product of $n-1$ matrices $\tilde{M}$
multiplying the initial length-1 vector with entry $f_{1}(0)=1$. Transferring
to the $\Omega$-notation by using Eqns. (24) and (25), we obtain

$\boxed{K_{n}\left(\Omega_{1},...,\Omega_{n-1}\right)=\\!\\!\\!\sum_{i_{n-1}=0}^{n-1}\sum_{i_{n-2}=0}^{n-2}\\!\cdots\\!\sum_{i_{2}=0}^{2}\sum_{i_{1}=0}^{1}\frac{\beta^{i_{n-1}}}{i_{n-1}+1}M_{1}(i_{n-1},i_{n-2})M_{2}(i_{n-2},i_{n-3})\cdots
M_{n-2}(i_{2},i_{1})M_{n-1}(i_{1},0)}$ (26)

with

$M_{j}(l,p)\equiv\frac{p!}{l!}\left[\delta_{l,p+1}\delta_{j}-\theta\left(p-l+1/2\right)\Delta_{j}^{1+p-l}\right].$
(27)

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).

##  V Explicit kernel functions
$K_{n}\bigl{(}\Omega_{1},...,\Omega_{n-1}\bigr{)}$ 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}\left(\Omega_{1},...,\Omega_{n-1}\right)$ 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 denoted 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 $\delta_{j}$ and
$\Delta_{j}$ for $j=1,2,...,n-1$ which appear in Tab. 1 are defined by

$\delta_{j}$ $\displaystyle\equiv$
$\displaystyle\delta_{\Omega_{1}+\Omega_{2}+...+\Omega_{j},0},$ (28)
$\Delta_{j}$ $\displaystyle\equiv$
$\displaystyle\Delta_{\Omega_{1}+\Omega_{2}+...+\Omega_{j}}\equiv\begin{cases}0&\textnormal{{if}
$\Omega_{1}$+$\Omega_{2}$+...+$\Omega_{j}$=0}\\\
\frac{1}{\Omega_{1}+\Omega_{2}+...+\Omega_{j}}&\textnormal{{if}
$\Omega_{1}$+$\Omega_{2}$+...+$\Omega_{j}$$\neq 0$}\end{cases},$ (29)

compare also to the previous section. As an example, for $n=2$ and $n=3$ we
obtain from Tab. 1

$\displaystyle K_{2}(\Omega_{1})$ $\displaystyle=$
$\displaystyle-\Delta_{1}+\frac{\beta}{2}\delta_{1},$ (30) $\displaystyle
K_{3}(\Omega_{1},\Omega_{2})$ $\displaystyle=$
$\displaystyle+\Delta_{1}\Delta_{2}-\frac{\beta}{2}\delta_{1}\Delta_{2}-\Delta_{1}\delta_{2}\left(\frac{\beta}{2}+\Delta_{1}\right)+\delta_{1}\delta_{2}\frac{\beta}{2}\frac{\beta}{3},$
(31)

respectively. The rows of Tab. 1 are organized with respect to the number $a$
of factors $\delta_{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\leq 3$
(Shvaika, 2006) or to fermionic $n=4$ (Hafermann et al., 2009; Shvaika, 2016;
Kugler et al., 2021) with $a=0,1$ (and $a=2,3$ guaranteed to vanish, see
below) or arbitrary $n$ with $a=0$ (Kugler et al., 2021). Alternative
expressions for the $n=3,4$ kernel functions with $a\leq 1$ were given in
(Kugler et al., 2021), but they are consistent with our kernel functions as
they yield the same correlation functions, see the Appendix.

#anom. | $K_{2}(\Omega_{1})$ | $K_{3}(\Omega_{1},\Omega_{2})$ | $K_{4}(\Omega_{1},\Omega_{2},\Omega_{3})$ | $K_{5}(\Omega_{1},\Omega_{2},\Omega_{3},\Omega_{4})$ | factor for entire row
---|---|---|---|---|---
$a=0$ | $-\Delta_{1}$ | $+\Delta_{1}\Delta_{2}$ | $-\Delta_{1}\Delta_{2}\Delta_{3}$ | $+\Delta_{1}\Delta_{2}\Delta_{3}\Delta_{4}$ | $1$
$a=1$ | $+\delta_{1}$ | $-\delta_{1}\Delta_{2}$ | $+\delta_{1}\Delta_{2}\Delta_{3}$ | $-\delta_{1}\Delta_{2}\Delta_{3}\Delta_{4}$ | $\frac{\beta}{2}$
| $-\Delta_{1}\delta_{2}$ | $+\Delta_{1}\delta_{2}\Delta_{3}$ | $-\Delta_{1}\delta_{2}\Delta_{3}\Delta_{4}$ | $\frac{\beta}{2}+\Delta_{1}$
|  | $+\Delta_{1}\Delta_{2}\delta_{3}$ | $-\Delta_{1}\Delta_{2}\delta_{3}\Delta_{4}$ | $\frac{\beta}{2}+\Delta_{1}+\Delta_{2}$
|  |  | $-\Delta_{1}\Delta_{2}\Delta_{3}\delta_{4}$ | $\frac{\beta}{2}+\Delta_{1}+\Delta_{2}+\Delta_{3}$
$a=2$ |  | $+\delta_{1}\delta_{2}$ | $-\delta_{1}\delta_{2}\Delta_{3}$ | $+\delta_{1}\delta_{2}\Delta_{3}\Delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}$
|  | $-\delta_{1}\Delta_{2}\delta_{3}$ | $+\delta_{1}\Delta_{2}\delta_{3}\Delta_{4}$ | $\frac{\beta}{2}\left(\frac{\beta}{3}+\Delta_{2}\right)$
|  | $-\Delta_{1}\delta_{2}\delta_{3}$ | $+\Delta_{1}\delta_{2}\delta_{3}\Delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}+\Delta_{1}\left(\frac{\beta}{2}+\Delta_{1}\right)$
|  |  | $+\delta_{1}\Delta_{2}\Delta_{3}\delta_{4}$ | $\frac{\beta}{2}\left(\frac{\beta}{3}+\Delta_{2}+\Delta_{3}\right)$
|  |  | $+\Delta_{1}\delta_{2}\Delta_{3}\delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}+\left(\Delta_{1}+\Delta_{3}\right)\left(\frac{\beta}{2}+\Delta_{1}\right)$
|  |  | $+\Delta_{1}\Delta_{2}\delta_{3}\delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}+\left(\Delta_{1}+\Delta_{2}\right)\left(\frac{\beta}{2}+\Delta_{2}\right)+\Delta_{1}^{2}$
$a=3$ |  |  | $+\delta_{1}\delta_{2}\delta_{3}$ | $-\delta_{1}\delta_{2}\delta_{3}\Delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}\frac{\beta}{4}$
|  |  | $-\delta_{1}\delta_{2}\Delta_{3}\delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}\left(\frac{\beta}{4}+\Delta_{3}\right)$
|  |  | $-\delta_{1}\Delta_{2}\delta_{3}\delta_{4}$ | $\frac{\beta}{2}\left(\frac{\beta}{3}\frac{\beta}{4}+\Delta_{2}\left(\frac{\beta}{3}+\Delta_{2}\right)\right)$
|  |  | $-\Delta_{1}\delta_{2}\delta_{3}\delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}\frac{\beta}{4}+\Delta_{1}\left(\frac{\beta}{2}\frac{\beta}{3}+\Delta_{1}\left(\frac{\beta}{2}+\Delta_{1}\right)\right)$
$a=4$ |  |  |  | $+\delta_{1}\delta_{2}\delta_{3}\delta_{4}$ | $\frac{\beta}{2}\frac{\beta}{3}\frac{\beta}{4}\frac{\beta}{5}$
Table 1: Universal kernel functions
$K_{n}\bigl{(}\Omega_{1},...,\Omega_{n-1}\bigr{)}$ for $n=2,3,4,5$ defined in
Eq. (14) and calculated from Eq. (26) in Sec. IV. 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
$\delta_{j}$ and $\Delta_{j}$ are defined in Eqns. (28) and (29). The rows are
organized with respect to the number $a$ of appearances of $\delta_{j}$, i.e.
the order of the anomalous terms.

In the case of purely fermionic correlators (all $A_{k}$ fermionic),
individual Matsubara frequencies $\omega_{k}$ cannot be zero. Thus the complex
numbers
$\Omega_{k}^{\underline{a}\,\underline{b}}=i\omega_{k}+E_{\underline{a}}-E_{\underline{b}}$
of Eq. (12) always have a finite imaginary part, regardless of the
eigenenergies. In this case, only sums of an even number of frequencies can be
zero, and we can simplify $\delta_{1}=\delta_{3}=\delta_{5}=...=0$. The
expressions for the kernels in Tab. 1, now denoted by $K_{n}|_{F}$ for the
fermionic case, simplify to

$\displaystyle K_{2}(\Omega_{1})|_{F}$ $\displaystyle=$
$\displaystyle-\Delta_{1},$ (32) $\displaystyle
K_{4}(\Omega_{1},\Omega_{2},\Omega_{3})|_{F}$ $\displaystyle=$
$\displaystyle\Delta_{1}\Delta_{3}\left[\delta_{2}\left(\frac{\beta}{2}+\Delta_{1}\right)-\Delta_{2}\right],$
(33) $\displaystyle K_{6}(\Omega_{1},...,\Omega_{5})|_{F}$ $\displaystyle=$
$\displaystyle\Delta_{1}\Delta_{3}\Delta_{5}\biggl{\\{}-\Delta_{2}\Delta_{4}-\delta_{2}\delta_{4}\left[\frac{\beta}{2}\frac{\beta}{3}+\left(\Delta_{1}+\Delta_{3}\right)\left(\frac{\beta}{2}+\Delta_{1}\right)\right]$
$\displaystyle+\delta_{4}\Delta_{2}\left(\frac{\beta}{2}+\Delta_{1}+\Delta_{2}+\Delta_{3}\right)+\delta_{2}\Delta_{4}\left(\frac{\beta}{2}+\Delta_{1}\right)\biggr{\\}}.$

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.

##  VI Applications: Hubbard atom and Free Spin

###  VI.1 Fermionic Hubbard atom

The Hubbard atom (HA) describes an isolated impurity or otherwise localized
system with Hamiltonian

$H=\epsilon(n_{\uparrow}+n_{\downarrow})+Un_{\uparrow}n_{\downarrow}-h(n_{\uparrow}-n_{\downarrow}),$
(35)

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_{\sigma}=d_{\sigma}^{\dagger}d_{\sigma}$ count the number of fermionic
particles with spin $\sigma\in\\{\uparrow,\downarrow\\}$, each contributing an
onsite energy $\epsilon$ 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 (Krien and Valli,
2019; Krien et al., 2021; Kappl et al., 2023). 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 divergences of two-
line irreducible vertex functions (Schäfer et al., 2016; Thunström et al.,
2018; Chalupa et al., 2018; Pelz et al., 2023) and signatures of the local
moment formation in generalized susceptibilities (Chalupa et al., 2021; Adler
et al., 2022). Using the fermionic kernels in Eqns. (32) and (33), we have
checked that our formalism reproduces the results for the 2-point and 4-point
correlators given in Refs. (Hafermann et al., 2009; Rohringer, 2013; Kugler et
al., 2021) for half-filling, $\epsilon=-U/2$ and $h=0$.

Correlation functions including bosonic operators describe the asymptotic
behaviour of the $n=4$ fermion vertex for large frequencies (Wentzell et al.,
2020) or the interaction of electrons by the exchange of an effective boson
(Krien et al., 2019; Gievers et al., 2022). 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\leq 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 (Kappl et
al., 2023). 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 $\delta_{\epsilon}n$ or take the form of a magnetic field
$\boldsymbol{\delta}_{h}\cdot\mathbf{S}$. Here, the parameters
$\delta_{\varepsilon}$ and $\boldsymbol{\delta}_{h}$ denote the strength of
the perturbation and we define

$\begin{aligned} n&=n_{\uparrow}+n_{\downarrow}\end{aligned},\quad
S^{x}=\frac{1}{2}\left(d_{\uparrow}^{\dagger}d_{\downarrow}+d_{\downarrow}^{\dagger}d_{\uparrow}\right),\quad
S^{y}=\frac{-i}{2}\left(d_{\uparrow}^{\dagger}d_{\downarrow}-d_{\downarrow}^{\dagger}d_{\uparrow}\right),\quad
S^{z}=\frac{1}{2}\left(n_{\uparrow}-n_{\downarrow}\right).$ (36)

The resulting changes of the expectation values of the density or
magnetization in arbitrary direction are described in second order of the
perturbation by the connected parts of the correlation functions
$G_{A_{1}A_{2}A_{3}}(\tau_{1},\tau_{2},\tau_{3})$, with
$A_{i}\in\\{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. (Kappl et al., 2023). In the following, we give explicit,
analytic expressions of the full correlation functions
$G_{A_{1}A_{2}A_{3}}(\omega_{1},\omega_{2})$ (i.e. including disconnected
parts), for arbitrary parameters $\epsilon$, $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. 1(b)] describe an empty
($|0\rangle$), singly occupied
($d_{\uparrow}^{\dagger}|0\rangle=|\uparrow\rangle$,
$d_{\downarrow}^{\dagger}|0\rangle=|\downarrow\rangle$) or doubly occupied
($d_{\uparrow}^{\dagger}d_{\downarrow}^{\dagger}|0\rangle=|\uparrow\downarrow\rangle$)
impurity with eigenenergies $E_{0}=0$, $E_{\uparrow}=\epsilon-h$,
$E_{\downarrow}=\epsilon+h$ and $E_{\uparrow\downarrow}=2\epsilon+U$,
respectively. The partition function is
$Z=1+e^{-\beta(\epsilon-h)}+e^{-\beta(\epsilon+h)}+e^{-\beta(2\epsilon+U)}.$
We define

$s=\frac{e^{-\beta\epsilon}}{Z}\sinh(\beta h),\quad
c=\frac{e^{-\beta\epsilon}}{Z}\cosh(\beta h),$ (37)

and obtain all non-vanishing bosonic 3-point correlation functions (where
$\omega_{3}=-\omega_{1}-\omega_{2}$):

$\displaystyle G_{nnn}(\omega_{1},\omega_{2})$
$\displaystyle=2\beta^{2}\delta_{\omega_{1}}\delta_{\omega_{2}}\left(\frac{4e^{-\beta(2\epsilon+U)}}{Z}+c\right),$
(38) $\displaystyle G_{nnS^{z}}(\omega_{1},\omega_{2})$
$\displaystyle=\beta^{2}\delta_{\omega_{1}}\delta_{\omega_{2}}s,$ (39)
$\displaystyle G_{nS^{x}S^{y}}(\omega_{1},\omega_{2})$
$\displaystyle=-\beta\delta_{\omega_{1}}s\frac{\omega_{2}}{\omega_{2}^{2}+4h^{2}},$
(40) $\displaystyle
G_{nS^{x}S^{x}}(\omega_{1},\omega_{2})=G_{nS^{y}S^{y}}(\omega_{1},\omega_{2})$
$\displaystyle=2\beta\delta_{\omega_{1}}\frac{h\ s}{\omega_{2}^{2}+4h^{2}},$
(41) $\displaystyle G_{nS^{z}S^{z}}(\omega_{1},\omega_{2})$
$\displaystyle=\frac{\beta^{2}}{2}\delta_{\omega_{1}}\delta_{\omega_{2}}c,$
(42) $\displaystyle
G_{S^{z}S^{x}S^{x}}(\omega_{1},\omega_{2})=G_{S^{z}S^{y}S^{y}}(\omega_{1},\omega_{2})$
$\displaystyle=-s\frac{\omega_{2}\omega_{3}+4h^{2}}{(\omega_{2}^{2}+4h^{2})(\omega_{3}^{2}+4h^{2})}+\beta\delta_{\omega_{1}}\frac{h\
c}{\omega_{2}^{2}+4h^{2}},$ (43) $\displaystyle
G_{S^{z}S^{z}S^{z}}(\omega_{1},\omega_{2})$
$\displaystyle=\frac{\beta^{2}}{4}\delta_{\omega_{1}}\delta_{\omega_{2}}s,$
(44) $\displaystyle G_{S^{x}S^{y}S^{z}}(\omega_{1},\omega_{2})$
$\displaystyle=2h\
s\frac{\omega_{1}-\omega_{2}}{(\omega_{1}^{2}+4h^{2})(\omega_{2}^{2}+4h^{2})}-\frac{\beta}{2}\delta_{\omega_{3}}c\frac{\omega_{1}}{\omega_{1}^{2}+4h^{2}}.$
(45)

We observe that each conserved quantity, in this case $n$ and $S_{z}$,
contributes an anomalous term $\propto\delta_{\omega_{k}}$in its respective
frequency argument $\omega_{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}^{\underline{1}\,\underline{2}}=A_{k}^{\underline{1}\,\underline{1}}\delta_{\underline{1},\underline{2}}$.
If $A_{k}^{\underline{1}\,\underline{1}}\neq 0$ for some state $\underline{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\rightarrow 0$, we introduce an additional
degeneracy $E_{\uparrow}=E_{\downarrow}=\epsilon$ 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

$\lim_{h\rightarrow 0}\frac{h\ \sinh(\beta
h)}{\omega_{k}^{2}+4h^{2}}=\frac{\beta}{4}\delta_{\omega_{k}},$ (46)

resulting in

$\displaystyle G_{nnn}(\omega_{1},\omega_{2})$
$\displaystyle=\beta^{2}\delta_{\omega_{1}}\delta_{\omega_{2}}\frac{2(4e^{-\beta(2\epsilon+U)}+e^{-\beta\epsilon})}{Z},$
(47) $\displaystyle G_{nS^{\alpha}S^{\alpha}}(\omega_{1},\omega_{2})$
$\displaystyle=\beta^{2}\delta_{\omega_{1}}\delta_{\omega_{2}}\frac{e^{-\beta\epsilon}}{2Z}\quad(\alpha\in\\{x,y,z\\}),$
(48) $\displaystyle G_{S^{x}S^{y}S^{z}}(\omega_{1},\omega_{2})$
$\displaystyle=\beta\frac{e^{-\beta\epsilon}}{2Z}(-\delta_{\omega_{1}}\Delta_{\omega_{2}}+\delta_{\omega_{2}}\Delta_{\omega_{1}}-\delta_{\omega_{1}+\omega_{2}}\Delta_{\omega_{1}}),$
(49)

with all other correlation functions vanishing. As already pointed out in Ref.
(Kappl et al., 2023), only the last correlation function retains a nontrivial
frequency dependence due to non-commuting operators.

###  VI.2 Free spin $S$

We now consider correlation functions of a free spin of length $S$, without a
magnetic field, so that temperature $T=1/\beta$ is the only energy scale. The
operators $\\{S^{\alpha}\\}_{\alpha=x,y,z}$ fulfill
$S^{x}S^{x}+S^{y}S^{y}+S^{z}S^{z}=S(S+1)$ and the SU(2) algebra
$\left[S^{\alpha_{1}},S^{\alpha_{2}}\right]=i\sum_{\alpha_{3}=\\{x,y,z\\}}\epsilon^{\alpha_{1}\alpha_{2}\alpha_{3}}S^{\alpha_{3}}$,
thus they are bosonic. Since the Hamiltonian vanishes and therefore all
eigenenergies are zero, every $\Omega_{k}^{\underline{a}\,\underline{b}}$ 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^{\alpha}(\tau)=S^{\alpha}$, the non-trivial frequency dependence of the
correlators, which can be can be non-vanishing at any order $n$, 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. (Krieg
and Kopietz, 2019; Goll et al., 2019, 2020; Tarasevych and Kopietz, 2021;
Tarasevych et al., 2022). 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^{\alpha_{1}}=S^{\alpha_{2}}=...=S^{\alpha_{n}}$
where the SU(2) algebra does not matter, or for finite magnetic field via an
equation of motion (Goll et al., 2019) or diagrammatic approach (Vaks and
Pikin, 1968; Vaks et al., 1968).

We define the spin raising and lowering operators,

$S^{\pm}=\left(S^{x}\pm iS^{y}\right)/\sqrt{2},$ (50)

which have to appear in pairs for a non-vanishing correlator due to spin-
rotation symmetry. As for the HA, we do not consider connected correlators in
this work for brevity. The classical $S^{z}$-correlator can be found from its
generating functional with source field $h$ (Krieg and Kopietz, 2019),

$\displaystyle\mathcal{G}\left(y=\beta h\right)$ $\displaystyle=$
$\displaystyle\frac{\sinh\left[y(S+1/2)\right]}{(2S+1)\sinh\left[y/2\right]},$
(51) $\displaystyle\left\langle(S^{z})^{l}\right\rangle$ $\displaystyle=$
$\displaystyle\underset{y\rightarrow
0}{\mathrm{lim}}\partial_{y}^{l}\mathcal{G}(y)\equiv b_{l-1},$ (52)

for example $b_{1}=\frac{S}{3}(S+1)$ and
$b_{3}=\frac{S}{15}\left(3S^{3}+6S^{2}+2S-1\right)$ and vanishing $b_{l}$ for
even $l$. For all other correlators involving $\alpha_{k}=\pm$, we adapt Eq.
(15) for the free spin case,

$G_{S^{\alpha_{1}}S^{\alpha_{2}}...S^{\alpha_{n}}}\left(\omega_{1},...,\omega_{n-1}\right)=\sum_{p\in
S_{n}}\left\langle
S^{\alpha_{p(1)}}S^{\alpha_{p(2)}}...S^{\alpha_{p(n)}}\right\rangle
K_{n}\left(i\omega_{p(1)},i\omega_{p(2)},...,i\omega_{p(n-1)}\right),$ (53)

where we made use of the fact that all eigenenergies are zero and the
Heisenberg time evolution is trivial. It is convenient to evaluate the equal-
time correlators in Eq. (53) as

$\left\langle
S^{\alpha_{1}}S^{\alpha_{2}}...S^{\alpha_{n}}\right\rangle=\frac{1}{2S+1}\sum_{m=-S}^{S}\left\langle
m\right|S^{\alpha_{1}}S^{\alpha_{2}}...S^{\alpha_{n}}\left|m\right\rangle\equiv\frac{1}{2S+1}\sum_{m=-S}^{S}\sum_{l=0}^{n}p_{l}m^{l}=p_{0}+\sum_{l=2}^{n}p_{l}b_{l-1}$
(54)

where in the last step we used Eq. (52). We find the real expansion
coefficients $\\{p_{l}\\}_{l=0,1,...,n}$ iteratively by moving through the
string $\alpha_{1}\alpha_{2}...\alpha_{n}$ from the right and start from
$p_{l}=\delta_{0,l}$. Based on the $S^{z}$ eigenstates
$\\{\bigl{|}m\bigr{\rangle}\\}_{m=-S,...,S-1,S}$ we obtain the iteration rules
from $S^{z}\bigl{|}m\bigr{\rangle}=m\bigl{|}m\bigr{\rangle}$ and
$S^{\pm}\bigl{|}m\bigr{\rangle}=\sqrt{1/2}\sqrt{S(S+1)-m(m\pm 1)}\bigl{|}m\pm
1\bigr{\rangle}$. We define an auxiliary integer $c$ that keeps track of the
intermediate state $\left|m+c\right\rangle$, initially $c=0$. Depending on the
$\alpha_{j}$ that we find in step $j=n,n-1...,1$ we take one of the following
actions: (i) For $\alpha_{j}=z$, we update $p_{l}\leftarrow
p_{l-1}+cp_{l}\;\forall l$ and leave $c$ unchanged. It is understood that
$p_{l<0}=0$. (ii) For $\alpha_{j}=+$, we combine the square-root factor
brought by the raising operator with the factor that comes from the necessary
$\alpha_{j^{\prime}}=-$ at another place in the string. We replace
$p_{l}\leftarrow-\frac{1}{2}p_{l-2}-\frac{2c+1}{2}p_{l-1}+\left(\frac{3}{2}b_{1}-c\frac{c+1}{2}\right)p_{l}\;\forall
l$ and then let $c\leftarrow c+1$. (iii) For $\alpha_{j}=-$, we update
$c\leftarrow c-1$ and keep $p_{l}$ unchanged, $p_{l}\leftarrow p_{l}\;\forall
l$.

$n=2$ | $G_{S^{+}S^{-}}(\omega)=G_{S^{z}S^{z}}(\omega)=\beta\delta_{\omega}b_{1}$
---|---
$n=3$ | $G_{S^{+}S^{-}S^{z}}(\omega_{1},\omega_{2})=\beta b_{1}(-\delta_{\omega_{1}}\Delta_{i\omega_{2}}+\delta_{\omega_{2}}\Delta_{i\omega_{1}}+\delta_{\omega_{1}+\omega_{2}}\Delta_{i\omega_{2}})=-iG_{S^{x}S^{y}S^{z}}(\omega_{1},\omega_{2})$
$n=4$ | $G_{S^{z}S^{z}S^{z}S^{z}}\left(\omega_{1},\omega_{2},\omega_{3}\right)=\delta_{\omega_{1}}\delta_{\omega_{2}}\delta_{\omega_{3}}\beta^{3}b_{3}$
$G_{S^{+}S^{+}S^{-}S^{-}}\left(\omega_{1},\omega_{2},\omega_{3}\right)=\beta
b_{1}[2\times\delta_{\omega_{1}}\delta_{\omega_{2}}\delta_{\omega_{3}}\times\frac{\beta^{2}}{5}\left(3b_{1}-\frac{1}{3}\right)+r]$
$G_{S^{+}S^{-}S^{z}S^{z}}\,\left(\omega_{1},\omega_{2},\omega_{3}\right)\,=\beta
b_{1}[1\times\delta_{\omega_{1}}\delta_{\omega_{2}}\delta_{\omega_{3}}\times\frac{\beta^{2}}{5}\left(3b_{1}-\frac{1}{3}\right)-r]$
$r=\Delta_{i\omega_{1}}\Delta_{i\text{$\omega_{2}$}}\left(\delta_{\omega_{1}+\omega_{3}}+\delta_{\omega_{2}+\omega_{3}}-\delta_{\omega_{3}}-\delta_{\omega_{4}}\right)-\left(\delta_{\omega_{1}}\Delta_{i\text{$\omega_{2}$}}^{2}+\delta_{\omega_{2}}\Delta_{i\text{$\omega_{1}$}}^{2}\right)\left(\delta_{\omega_{3}}+\delta_{\text{$\omega_{4}$}}\right)-\Delta_{i\omega_{3}}\Delta_{i\omega_{4}}\left(\delta_{\omega_{1}}+\delta_{\omega_{2}}\right)$
Table 2: Matsubara correlation functions for a free spin-$S$ up to order
$n=4$. Here, $\omega_{4}=-\omega_{1}-\omega_{2}-\omega_{3}$.

Our final results for the free spin correlators are reported in Tab. 2. We
reproduce the known spin correlators for $n=2,3$ and determine the non-
classical correlators $G_{S^{+}S^{+}S^{-}S^{-}}$ and
$G_{S^{+}S^{-}S^{z}S^{z}}$ at order $n=4$, which to the best of our knowledge
were not available in the literature (tha, ). 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

$\displaystyle\Delta_{a+b}\left(\Delta_{a}+\Delta_{b}\right)-\Delta_{a}\Delta_{b}$
$\displaystyle=$
$\displaystyle\delta_{a}\Delta_{b}^{2}+\delta_{b}\Delta_{a}^{2}-\delta_{a+b}\Delta_{a}\Delta_{b}.$
(55)

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
$\bigl{|}\uparrow\bigr{\rangle},\bigl{|}\downarrow\bigr{\rangle}$. 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^{-\beta\epsilon}/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 local-moment regime.

##  VII 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-of-motion 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\rightarrow 0$ where the
formally divergent anomalous contributions are to be understood as
$\beta\delta_{\omega,0}\rightarrow 2\pi\delta(\omega)$. 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
(Rohringer, 2013). 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 (Ge et al., ).

###### Acknowledgements.

We acknowledge useful discussions with Karsten Held, Friedrich Krien, Seung-
Sup Lee, Peter Kopietz, Fabian Kugler, Nepomuk Ritz, Georg Rohringer, Andreas
Rückriegel. We thank Andreas Rückriegel for sharing unpublished results on
4-point free spin correlators and pointing out further simplifications of the
analytical expressions. BS and BS are supported by a MCQST-START fellowship.
We acknowledge funding from the International Max Planck Research School for
Quantum Science and Technology (IMPRSQST) for JH, from the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy EXC-2111 (Project
No. 390814868), and from the Munich Quantum Valley, supported by the Bavarian
state government with funds from the Hightech Agenda Bayern Plus.

## Appendix: Equivalence to convention of Ref. [21]

In Ref. (Kugler et al., 2021) by Kugler, Lee and von Delft (KLD), only regular
($a=0$) and anomalous terms of order $a=1$ have been considered for $n=3,4$.
The corresponding kernel functions were derived from only $(n-1)!$
permutations by setting $\tau_{n}=0$ and $\tau_{i\neq n}>0$, but still applied
to all $n!$ permutations to obtain the correlation functions. For $n=3$, the
resulting kernel function (Eq. (46) in Ref. (Kugler et al., 2021)) reads

$K_{3,\textnormal{KLD}}(\Omega_{1},\Omega_{2})=\Delta_{1}\Delta_{2}-\Delta_{1}\delta_{2}\frac{1}{2}\left(\beta+\Delta_{1}\right)-\delta_{1}\Delta_{2}\frac{1}{2}\left(\beta+\Delta_{2}\right).$
(56)

This can be compared to the corresponding kernel function for $n=3$ found in
our Eq. (31) truncated to $a\leq 1$,

$K_{3}^{a\leq
1}(\Omega_{1},\Omega_{2})=\Delta_{1}\Delta_{2}-\Delta_{1}\delta_{2}\left(\frac{\beta}{2}+\Delta_{1}\right)-\frac{\beta}{2}\delta_{1}\Delta_{2}.$
(57)

Both approaches are equally valid and should yield the same correlation
functions (consistently discarding terms with $a=2$), yet the kernel functions
are obviously different. To resolve this issue, we define the difference of
the kernel functions

$K_{3,\textnormal{diff}}(\Omega_{1},\Omega_{2})=K_{3,\textnormal{KLD}}(\Omega_{1},\Omega_{2})-K_{3}^{a\leq
1}(\Omega_{1},\Omega_{2})=\frac{1}{2}\left(\Delta_{1}^{2}\delta_{2}-\delta_{1}\Delta_{2}^{2}\right)$
(58)

and show that the corresponding contributions to the correlation function
vanishes when summed over cyclically related permutations $p=123,231,312$.
These contributions are given by

$\displaystyle\frac{1}{Z}\sum_{p=123,231,312}\zeta(p)\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{p(1)}^{\underline{1}\underline{2}}A_{p(2)}^{\underline{2}\underline{3}}A_{p(3)}^{\underline{3}\underline{1}}K_{3,\textnormal{diff}}(\Omega_{p(1)}^{\underline{1}\,\underline{2}},\Omega_{p(2)}^{\underline{2}\,\underline{3}})$
(59)
$\displaystyle=\frac{\zeta(123)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{1}^{\underline{1}\underline{2}}A_{2}^{\underline{2}\underline{3}}A_{3}^{\underline{3}\underline{1}}\left(\frac{(1-\delta_{\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{2}}})\delta_{\omega_{1}+\omega_{2}}\delta_{E_{\underline{1}}-E_{\underline{3}}}}{(i\omega_{1}+E_{\underline{1}}-E_{\underline{2}})^{2}}-\frac{\delta_{\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{2}}}(1-\delta_{\omega_{1}+\omega_{2}}\delta_{E_{\underline{1}}-E_{\underline{3}}})}{(i\omega_{1}+i\omega_{2}+E_{\underline{1}}-E_{\underline{3}})^{2}}\right)$
$\displaystyle+\frac{\zeta(231)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{2}^{\underline{1}\underline{2}}A_{3}^{\underline{2}\underline{3}}A_{1}^{\underline{3}\underline{1}}\left(\frac{(1-\delta_{\omega_{2}}\delta_{E_{\underline{1}}-E_{\underline{2}}})\delta_{\omega_{2}+\omega_{3}}\delta_{E_{\underline{1}}-E_{\underline{3}}}}{(i\omega_{2}+E_{\underline{1}}-E_{\underline{2}})^{2}}-\frac{\delta_{\omega_{2}}\delta_{E_{\underline{1}}-E_{\underline{2}}}(1-\delta_{\omega_{2}+\omega_{3}}\delta_{E_{\underline{1}}-E_{\underline{3}}})}{(i\omega_{2}+i\omega_{3}+E_{\underline{1}}-E_{\underline{3}})^{2}}\right)$
$\displaystyle+\frac{\zeta(312)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{3}^{\underline{1}\underline{2}}A_{1}^{\underline{2}\underline{3}}A_{2}^{\underline{3}\underline{1}}\left(\frac{(1-\delta_{\omega_{3}}\delta_{E_{\underline{1}}-E_{\underline{2}}})\delta_{\omega_{3}+\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{3}}}}{(i\omega_{3}+E_{\underline{1}}-E_{\underline{2}})^{2}}-\frac{\delta_{\omega_{3}}\delta_{E_{\underline{1}}-E_{\underline{2}}}(1-\delta_{\omega_{3}+\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{3}}})}{(i\omega_{3}+i\omega_{1}+E_{\underline{1}}-E_{\underline{3}})^{2}}\right).$

Considering the second term of permutation $p=312$ and renaming the summation
variables $\underline{2}\rightarrow\underline{1}$,
$\underline{3}\rightarrow\underline{2}$,
$\underline{1}\rightarrow\underline{3}$ yields

$\displaystyle-\frac{\zeta(312)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{3}^{\underline{1}\underline{2}}A_{1}^{\underline{2}\underline{3}}A_{2}^{\underline{3}\underline{1}}\frac{\delta_{\omega_{3}}\delta_{E_{\underline{1}}-E_{\underline{2}}}(1-\delta_{\omega_{3}+\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{3}}})}{i(\omega_{3}+\omega_{1})+E_{\underline{1}}-E_{\underline{3}}}$
(60)
$\displaystyle=-\frac{\zeta(312)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{3}}}A_{1}^{\underline{1}\underline{2}}A_{2}^{\underline{2}\underline{3}}A_{3}^{\underline{3}\underline{1}}\frac{\delta_{\omega_{3}}\delta_{E_{\underline{3}}-E_{\underline{1}}}(1-\delta_{\omega_{3}+\omega_{1}}\delta_{E_{\underline{3}}-E_{\underline{2}}})}{(i\omega_{3}+i\omega_{1}+E_{\underline{3}}-E_{\underline{2}})^{2}}$
$\displaystyle=-\frac{\zeta(123)}{2Z}\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{1}^{\underline{1}\underline{2}}A_{2}^{\underline{2}\underline{3}}A_{3}^{\underline{3}\underline{1}}\frac{\delta_{\omega_{1}+\omega_{2}}\delta_{E_{\underline{1}}-E_{\underline{3}}}(1-\delta_{\omega_{1}}\delta_{E_{\underline{1}}-E_{\underline{2}}})}{(i\omega_{1}+E_{\underline{1}}-E_{\underline{2}})^{2}},$

where we used $\omega_{3}=-\omega_{1}-\omega_{2}$ and the fact that
$\delta_{\omega_{3}}$ enforces the third operator to be bosonic, such that
$\zeta(312)=\zeta(123)$. This term exactly cancels the first contribution of
permutation $p=123$ in (59). Repeating similar steps for the remaining terms,
we find that 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,
leading to

$\frac{1}{Z}\sum_{p\in\\{123,231,312\\}}\zeta(p)\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{p(1)}^{\underline{1}\underline{2}}A_{p(2)}^{\underline{2}\underline{3}}A_{p(3)}^{\underline{3}\underline{1}}K_{3,\textnormal{diff}}(\Omega_{p(1)}^{\underline{1}\underline{2}},\Omega_{p(2)}^{\underline{2}\underline{3}})=0.$
(61)

Similarly, summing over the second set of cyclically related permutations
$p=132,213,321$ leads to a vanishing result, leading to the conclusion that

$\frac{1}{Z}\sum_{p\in
S_{3}}\zeta(p)\sum_{\underline{1}\underline{2}\underline{3}}e^{-\beta
E_{\underline{1}}}A_{p(1)}^{\underline{1}\underline{2}}A_{p(2)}^{\underline{2}\underline{3}}A_{p(3)}^{\underline{3}\underline{1}}K_{3,\textnormal{diff}}(\Omega_{p(1)}^{\underline{1}\underline{2}},\Omega_{p(2)}^{\underline{2}\underline{3}})=0.$
(62)

Thus we have shown that both kernel functions in Eqns. (56) and (57) 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$.

## References

  * Negele and Orland [1988] J. W. Negele and H. Orland.  _Quantum Many-Particle Systems_.  Westview, 1988.
  * Bruus and Flensberg [2004] H. Bruus and K. Flensberg.  _Many-Body Quantum Theory in Condensed Matter Physics_.  Oxford Graduate Texts, 2004.
  * Kappl et al. [2023] P. Kappl, F. Krien, C. Watzenböck, and K. Held.  Non-linear responses and three-particle correlators in correlated electron systems exemplified by the Anderson impurity model.  _Phys. Rev. B_ , 107(20):205108, May 2023.  doi: 10.1103/PhysRevB.107.205108.
  * Semeghini et al. [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin.  Probing topological spin liquids on a programmable quantum simulator.  _Science_ , 374(6572):1242–1247, December 2021\.  doi: 10.1126/science.abi8794.
  * Hedin [1965] L. Hedin.  New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem.  _Phys. Rev._ , 139(3A):A796–A823, August 1965\.  doi: 10.1103/PhysRev.139.A796.
  * [6] N.E. Bickers.  _Self-Consistent Many-Body Theory for Condensed Matter Systems_.  Springer.
  * Roulet et al. [1969] B. Roulet, J. Gavoret, and P. Nozières.  Singularities in the X-Ray Absorption and Emission of Metals. I. First-Order Parquet Calculation.  _Phys. Rev._ , 178(3):1072–1083, February 1969\.  doi: 10.1103/PhysRev.178.1072.
  * Kopietz et al. [2010] P. Kopietz, L. Bartosch, and F. Schütz.  _Introduction to the Functional Renormalization Group_.  Springer, 2010.
  * Metzner et al. [2012] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer.  Functional renormalization group approach to correlated fermion systems.  _Rev. Mod. Phys._ , 84:299, March 2012.  doi: 10.1103/RevModPhys.84.299.
  * Wetterich [1993] C. Wetterich.  Exact evolution equation for the effective potential.  _Physics Letters B_ , 301(1):90–94, February 1993\.  ISSN 0370-2693.  doi: 10.1016/0370-2693(93)90726-X.
  * Lee et al. [2021] S.-S. Lee, F. Kugler, and J. von Delft.  Computing Local Multipoint Correlators Using the Numerical Renormalization Group.  _Phys. Rev. X_ , 11(4):041007, October 2021.  doi: 10.1103/PhysRevX.11.041007.
  * Reuther and Thomale [2014] J. Reuther and R. Thomale.  Cluster functional renormalization group.  _Phys. Rev. B_ , 89(2):024412, January 2014.  doi: 10.1103/PhysRevB.89.024412.
  * Krieg and Kopietz [2019] J. Krieg and P. Kopietz.  Exact renormalization group for quantum spin systems.  _Phys. Rev. B_ , 99(6):060403, February 2019.  doi: 10.1103/PhysRevB.99.060403.
  * Rückriegel et al. [2022] A. Rückriegel, J. Arnold, R. Goll, and P. Kopietz.  Spin functional renormalization group for dimerized quantum spin systems.  _Phys. Rev. B_ , 105(22):224406, June 2022.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.105.224406.
  * Rohringer et al. [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held.  Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory.  _Rev. Mod. Phys._ , 90(2):025003, May 2018.  doi: 10.1103/RevModPhys.90.025003.
  * Chalupa et al. [2021] P. Chalupa, T. Schäfer, M. Reitner, D. Springer, S. Andergassen, and A. Toschi.  Fingerprints of the Local Moment Formation and its Kondo Screening in the Generalized Susceptibilities of Many-Electron Problems.  _Phys. Rev. Lett._ , 126(5):056403, February 2021\.  doi: 10.1103/PhysRevLett.126.056403.
  * Lehmann [1954] H. Lehmann.  Über Eigenschaften von Ausbreitungsfunktionen und Renormierungskonstanten quantisierter Felder.  _Nuovo Cim_ , 11(4):342–357, April 1954.  ISSN 1827-6121.  doi: 10.1007/BF02783624.
  * Shvaika [2006] A.M. Shvaika.  On the spectral relations for multitime correlation functions.  9(3(47)):447–458, 2006.
  * Hafermann et al. [2009] H. Hafermann, C. Jung, S. Brener, M. I. Katsnelson, A. N. Rubtsov, and A. I. Lichtenstein.  Superperturbation solver for quantum impurity models.  _Europhys. Lett._ , 85(2):27007, January 2009\.  ISSN 0295-5075, 1286-4854.  doi: 10.1209/0295-5075/85/27007.
  * Shvaika [2016] A.M. Shvaika.  Spectral properties of four-time fermionic Green’s functions.  19(3):33004, 2016.  doi: 10.5488/CMP.19.33004.
  * Kugler et al. [2021] F. B. Kugler, S.-S. B. Lee, and J. von Delft.  Multipoint Correlation Functions: Spectral Representation and Numerical Evaluation.  _Phys. Rev. X_ , 11(4):041006, October 2021.  doi: 10.1103/PhysRevX.11.041006.
  * Kwok and Schultz [1969] P. C. Kwok and T. D. Schultz.  Correlation functions and Green functions: Zero-frequency anomalies.  _J. Phys. C: Solid State Phys._ , 2(7):1196, July 1969.  ISSN 0022-3719.  doi: 10.1088/0022-3719/2/7/312.
  * Watzenböck et al. [2022] C. Watzenböck, M. Fellinger, K. Held, and A. Toschi.  Long-term memory magnetic correlations in the Hubbard model: A dynamical mean-field theory analysis.  _SciPost Phys._ , 12(6):184, June 2022.  ISSN 2542-4653.  doi: 10.21468/SciPostPhys.12.6.184.
  * Tsvelik [2007] A. Tsvelik.  _Quantum Field Theory in Condensed Matter Physics_.  Cambridge U. Press, 2007.
  * Rohringer et al. [2012] G. Rohringer, A. Valli, and A. Toschi.  Local electronic correlation at the two-particle level.  _Phys. Rev. B_ , 86(12):125114, September 2012\.  doi: 10.1103/PhysRevB.86.125114.
  * Rohringer [2013] G. Rohringer.  _New Routes towards a Theoretical Treatment of Nonlocal Electronic Correlations_.  Thesis, Technische Universität Wien, 2013.
  * Krien and Valli [2019] F. Krien and A. Valli.  Parquetlike equations for the Hedin three-leg vertex.  _Phys. Rev. B_ , 100(24):245147, December 2019\.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.100.245147.
  * Krien et al. [2021] F. Krien, A. Kauch, and K. Held.  Tiling with triangles: Parquet and G W $\gamma$ methods unified.  _Phys. Rev. Research_ , 3(1):013149, February 2021\.  ISSN 2643-1564.  doi: 10.1103/PhysRevResearch.3.013149.
  * Schäfer et al. [2016] T. Schäfer, S. Ciuchi, M. Wallerberger, P. Thunström, O. Gunnarsson, G. Sangiovanni, G. Rohringer, and A. Toschi.  Nonperturbative landscape of the Mott-Hubbard transition: Multiple divergence lines around the critical endpoint.  _Phys. Rev. B_ , 94(23):235108, December 2016\.  doi: 10.1103/PhysRevB.94.235108.
  * Thunström et al. [2018] P. Thunström, O. Gunnarsson, Sergio Ciuchi, and G. Rohringer.  Analytical investigation of singularities in two-particle irreducible vertex functions of the Hubbard atom.  _Phys. Rev. B_ , 98(23):235107, December 2018\.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.98.235107.
  * Chalupa et al. [2018] P. Chalupa, P. Gunacker, T. Schäfer, K. Held, and A. Toschi.  Divergences of the irreducible vertex functions in correlated metallic systems: Insights from the Anderson impurity model.  _Phys. Rev. B_ , 97(24):245136, June 2018.  doi: 10.1103/PhysRevB.97.245136.
  * Pelz et al. [2023] M. Pelz, S. Adler, M. Reitner, and A. Toschi.  The highly nonperturbative nature of the Mott metal-insulator transition: Two-particle vertex divergences in the coexistence region.  (arXiv:2303.01914), March 2023.  doi: 10.48550/arXiv.2303.01914.
  * Adler et al. [2022] S. Adler, F. Krien, P. Chalupa-Gantner, G. Sangiovanni, and A. Toschi.  Non-perturbative intertwining between spin and charge correlations: A "smoking gun" single-boson-exchange result.  (arXiv:2212.09693), December 2022.  doi: 10.48550/arXiv.2212.09693.
  * Wentzell et al. [2020] N. Wentzell, G. Li, A. Tagliavini, C. Taranto, G. Rohringer, K. Held, A. Toschi, and S. Andergassen.  High-frequency asymptotics of the vertex function: Diagrammatic parametrization and algorithmic implementation.  _Phys. Rev. B_ , 102(8):085106, August 2020.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.102.085106.
  * Krien et al. [2019] F. Krien, A. Valli, and M. Capone.  Single-boson exchange decomposition of the vertex function.  _Phys. Rev. B_ , 100(15):155149, October 2019\.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.100.155149.
  * Gievers et al. [2022] M. Gievers, E. Walter, A. Ge, J. von Delft, and F. B. Kugler.  Multiloop flow equations for single-boson exchange fRG.  _Eur. Phys. J. B_ , 95(7):108, July 2022.  ISSN 1434-6036.  doi: 10.1140/epjb/s10051-022-00353-6.
  * Goll et al. [2019] R. Goll, D. Tarasevych, J. Krieg, and P. Kopietz.  Spin functional renormalization group for quantum Heisenberg ferromagnets: Magnetization and magnon damping in two dimensions.  _Phys. Rev. B_ , 100(17):174424, November 2019\.  doi: 10.1103/PhysRevB.100.174424.
  * Goll et al. [2020] R. Goll, A. Rückriegel, and P. Kopietz.  Zero-magnon sound in quantum Heisenberg ferromagnets.  _Phys. Rev. B_ , 102(22):224437, December 2020\.  doi: 10.1103/PhysRevB.102.224437.
  * Tarasevych and Kopietz [2021] D. Tarasevych and P. Kopietz.  Dissipative spin dynamics in hot quantum paramagnets.  _Phys. Rev. B_ , 104(2):024423, July 2021.  doi: 10.1103/PhysRevB.104.024423.
  * Tarasevych et al. [2022] D. Tarasevych, A. Rückriegel, S. Keupert, V. Mitsiioannou, and P. Kopietz.  Spin functional renormalization group for the J1-J2-J3 quantum Heisenberg model.  _Phys. Rev. B_ , 106(17):174412, November 2022\.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.106.174412.
  * Vaks and Pikin [1968] V. G. Vaks and S. A. Pikin.  Thermodynamics of an ideal ferromagnetic substance.  _Soviet Journal of Experimental and Theoretical Physics_ , 26(1):188, 1968.
  * Vaks et al. [1968] V. G. Vaks, A. I. Larkin, and S. A. Pikin.  Spin Waves and Correlation Functions in a Ferromagnetic.  _Soviet Journal of Experimental and Theoretical Physics_ , 26:647, March 1968.  ISSN 1063-7761.
  * [43] We thank Andreas Rückriegel for sharing unpublished results on 4-point free spin correlators and pointing out further simplifications of the analytical expressions.
  * [44] A. Ge, J. Halbinger, S.-S. Lee, J. von Delft, and F. B. Kugler.  manuscript in preparation.