# Frustrated Quantum Spins at finite Temperature: Pseudo-Majorana functional
RG approach

Nils Niggemann  Dahlem Center for Complex Quantum Systems and Institut für
Theoretische Physik, Freie Universität Berlin, Arnimallee 14, 14195 Berlin,
Germany  Björn Sbierski  Department of Physics, University of California,
Berkeley, California 94720, USA  Johannes Reuther  Dahlem Center for Complex
Quantum Systems and Institut für Theoretische Physik, Freie Universität
Berlin, Arnimallee 14, 14195 Berlin, Germany  Helmholtz-Zentrum für
Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany

(August 28, 2024)

###### Abstract

The pseudofermion functional renormalization group (PFFRG) method has proven
to be a powerful numerical approach to treat frustrated quantum spin systems.
In its usual implementation, however, the complex fermionic representation of
spin operators introduces unphysical Hilbert space sectors which render an
application at finite temperatures inaccurate. In this work we formulate a
general functional renormalization group approach based on Majorana fermions
to overcome these difficulties. We, particularly, implement spin operators via
an $SO(3)$ symmetric Majorana representation which does not introduce any
unphysical states and, hence, remains applicable to quantum spin models at
finite temperatures. We apply this scheme, dubbed pseudo Majorana functional
renormalization group (PMFRG) method, to frustrated Heisenberg models on small
spin clusters as well as square and triangular lattices. Computing the finite
temperature behavior of spin correlations and thermodynamic quantities such as
free energy and heat capacity, we find good agreement with exact
diagonalization and the high-temperature series expansion down to moderate
temperatures. We observe a significantly enhanced accuracy of the PMFRG
compared to the PFFRG at finite temperatures. More generally, we conclude that
the development of functional renormalization group approaches with Majorana
fermions considerably extends the scope of applicability of such methods.

##  I Introduction

Finding numerical solutions of quantum many-body problems is one of the core
disciplines in modern condensed matter theory. In a wide range of physical
settings the problem amounts to analyze ground-state and finite-temperature
phases of a system of interacting spins on a lattice. Even though the
corresponding microscopic models are often conceptually simple, such as two-
body Heisenberg spin Hamiltonians, they may harbor a colorful range of
physical phenomena including exotic types of long-range orders [1], quantum
phase transitions [2, 3] or quantum spin liquids [4, 5, 6]. While quantum spin
phases are traditionally described in terms of broken or unbroken symmetries,
a more modern understanding also includes concepts such as long-range
entanglement or topological order [7] and reaches out to applications in the
context of quantum information processing [8].

Despite the shifts of focus which the field has gone through in the recent
decades, the accurate numerical treatment of interacting quantum spin systems
remains a highly challenging and longstanding problem. In fact, none of the
currently available numerical methods is able to ultimately determine the
eigenstates of a generic spin model. For example, quantum Monte Carlo methods
[9, 10] which enjoy the invaluable advantage that numerical errors are, in
principle, only of statistical nature, suffer from the infamous sign problem
when applied to frustrated spin systems. Similarly, density matrix
renormalization group, matrix product, and tensor network approaches [11, 12,
13, 14, 15] have made tremendous progress in recent years and are the
undisputed method of choice for a variety of spin systems (particularly in one
dimension). On the other hand, the scaling of the entanglement entropy poses a
serious challenge for such techniques in higher dimensions.

An alternative approach is based on functional renormalization group (FRG)
concepts [16, 17, 18] which are, in principle, oblivious to the system’s
dimensionality. In its standard fermionic formulation this technique has first
been applied in the context of electronic Hubbard-like models [19, 20, 21]
where it has become an established tool to describe competing types of long-
range orders. In addition, a more recently developed variant of the FRG [22]
specifically targets quantum spin systems. The key conceptual step of this
latter technique is to express the spin operators in terms of auxiliary
fermions [23], justifying the name pseudofermion functional renormalization
group (PFFRG). Within the last decade the PFFRG has been successfully applied
to a wide range of spin systems [22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63] and has constantly been extended
and generalized. Today, the PFFRG is, hence, remarkably flexible with a scope
of applicability comprising two dimensional [22, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 37, 38, 40, 41, 44, 45, 46, 48, 49, 51, 54, 55, 57, 59,
60, 61, 63] and three dimensional [36, 39, 42, 43, 47, 50, 52, 53, 55, 56, 58,
62] quantum spin systems on arbitrary lattices, including complex frustrated
and longer-range coupled networks [48, 49] with general isotropic or
anisotropic [54] two-body spin interactions. Further recent developments
concern the generalization to arbitrary spin magnitudes $S$ [41] or higher
spin symmetry groups $SU(N)$ [44, 45, 60] and, on a more technical level, the
implementation of multi-loop schemes [46, 62, 63].

Despite its success in accurately determining ground state spin correlations,
the PFFRG comes along with a well-known obstacle. The aforementioned
pseudofermionic description introduces an enlargement of the Hilbert space
associated with states that do not correspond to states of the physical spin
system. These unphysical states typically appear at energies above the ground
state energy of the spin system. Thus, on the level of zero-temperature
investigations, this problem has been argued to be rather mild and can be
treated by shifting unphysical states to higher energies [41]. In a recent
investigation of this problem, on the other hand, the average spin magnitude
within the PFFRG was found to differ from the theoretically expected result
even for higher loop orders [63]. More importantly, the enlarged Hilbert space
has so far prohibited an application to finite temperatures.

This work aims at resolving issues due to unphysical spin states by modifying
the PFFRG on a very fundamental level. Instead of using a complex fermionic
spin representation, we employ a certain, so-called $SO(3)$ Majorana fermion
rewriting of spin operators [64, 65] which does not generate unphysical states
but only introduces redundant Hilbert space sectors. This property
distinguishes it from other Majorana representations [75] and as such makes it
attractive as a first candidate for a Majorana-based spin FRG. We,
accordingly, dub our approach pseudo Majorana functional renormalization group
(PMFRG) method. This modification opens up various directions of
investigation: $(i)$ Most importantly, the PMFRG becomes applicable to finite
temperatures which only requires small methodological adjustments presented
below. $(ii)$ As a side product, we discuss how to calculate thermodynamic
quantities such as the free energy, energy and heat capacity which have so far
not been studied within the PFFRG. $(iii)$ To the best of our knowledge, a
Majorana-implementation of the FRG has so far not been published. Our
developments below are formulated in a general way such that they are
applicable to arbitrary Majorana models also outside the realm of quantum
magnetism. $(iv)$ Certain spin models, most prominently the Kitaev honeycomb
model [66], are exactly solvable when expressed in terms of Majorana fermions.
Although Kitaev’s spin representation differs from the one employed here, the
exact solution is also obtainable within the $SO(3)$ Majorana representation
[75] used here. Even though not the focus of this work, one may thus expect
that the PMFRG performs better for Kitaev-type spin models and perturbations
thereof as compared to the PFFRG.

Apart from the methodological focus of this work, we also present various
applications of the PMFRG to simple quantum spin models allowing us to assess
its accuracy. As a first benchmark test we treat small clusters of up to six
interacting spins where our results can be straightforwardly compared with
exact diagonalization. An overall finding is that the thermodynamic behavior
of the spin correlations from PMFRG are surprisingly accurate and reproduce
the exact result significantly better than PFFRG. It should be emphasized that
despite the finite Hilbert space of our spin clusters, their treatment within
PMFRG is still highly non-trivial and poses the same challenges as for
infinite lattice systems. Indeed, due to the incorporation of various mean-
field limits, one can expect that the FRG unfolds its full strength only in
infinite spin systems of two and higher dimensions. This motivates us to move
on to frustrated Heisenberg models on 2D square and triangular lattices where
we, likewise, find good agreement of thermodynamic properties with other
approaches. A persistent technical issue, however, occurs in the low
temperature limit where PMFRG detects spurious divergencies of spin
correlations. We interpret this behavior as an artifact of the redundant
Hilbert space sectors in our Majorana representation. While such subtleties
remain to be further studied we expect that our developments lay the
groundwork for various future directions of research and significantly enlarge
the scope of applicability of FRG approaches.

The remainder of this work is organized as follows: After briefly reviewing
the key concepts of the PFFRG in Sec. II, we discuss in detail the properties
of the $SO(3)$ Majorana representation in Sec. III. Thereafter, Sec. IV
formulates a general functional renormalization group approach for Majorana
systems. The specific implementation for Heisenberg spin models in $SO(3)$
Majorana representation is discussed in Section V with a particular focus on
the parametrization of vertex functions, taking into account the system’s
symmetries. The resulting RG flow equations are presented in Sec. VI and the
computation of various physical observables is detailed in Sec. VII. The
following Secs. VIII and IX discuss applications to small interacting spin
clusters as well as to square and triangular lattice models. The paper ends
with a conclusion in Sec. X.

##  II Basic concepts of the PFFRG

As a preparation for the following sections, we first briefly review basic
concepts and properties of the PFFRG approach without being exhaustive on all
methodological details. For a more detailed and self-contained description, we
refer the interested reader to Refs. [22, 41, 44, 54].

The PFFRG is capable of treating general two-body spin Hamiltonians; in this
work, however, only Heisenberg models of the form

$H=\sum_{(i,j)}J_{ij}\sum_{\alpha}S_{i}^{\alpha}S_{j}^{\alpha}$ (1)

will be considered, where $(i,j)$ refers to all possible pairings of sites and
$S^{\alpha}_{i}$ is the $\alpha$ component of a spin-$1/2$ operator at site
$i$. We note in passing that recently developed FRG approaches [67, 68]
directly take Eq. (1) as a starting point. In contrast, the PFFRG treats the
interacting fermionic model that results from representing the spin-1/2
operators via (pseudo-) fermions $f_{ia}$ (with $a=\uparrow,\downarrow$) [23]:

$S^{\alpha}_{i}=\frac{1}{2}\sum_{a,b}f_{ia}^{\dagger}\sigma^{\alpha}_{ab}f_{ib}\;.$
(2)

Here and in the following, we set $\hbar=k_{B}=1$. However, this
representation is a valid rewriting of the spin operators only in the local
subspace with $\sum_{a}f_{ia}^{\dagger}f_{ia}=1$ while states with zero or
double fermionic occupancy are unphysical. Since these spurious states carry
zero spin, they may be considered as voids in the spin system, associated with
an excitation energy on the order of the exchange coupling. As a consequence,
ground state properties are believed to be largely unaffected by unphysical
states, such that at $T=0$ the PFFRG may be faithfully implemented with the
simpler condition $\sum_{a}\langle f_{ia}^{\dagger}f_{ia}\rangle=1$. Other
approaches aiming to enforce the occupancy constraint more rigorously
introduce an energy penalty for unphysical states [41] or a particular form of
an imaginary chemical potential [69]. In either case, the unphysical states
remain an obstacle for an application of the PFFRG, especially at finite
temperatures. This motivates us to implement the FRG with the Majorana
representation discussed in Sec. III where no unphysical states occur.

The key benefit of the representation in Eq. (2) is that the resulting model
becomes amenable to fermionic many-body techniques such as the FRG which is
formulated in terms of irreducible fermionic vertex functions (“essential
parts of correlation functions”). The centerpiece of the method is given by a
hierarchy of flow equations reminiscent of one-loop diagrammatic perturbation
theory which describe the change of vertex functions when a Matsubara-
frequency cutoff parameter $\Lambda$, introduced in the bare Green function
$G^{0,\Lambda}(i\omega)=G^{0}(i\omega)\Theta(|\omega|-\Lambda)$, is varied.
The basic idea is that at the starting point $\Lambda=\infty$, the bare
propagator vanishes and all vertex functions are trivially known. For a
numerical solution of the flow equations down to $\Lambda=0$ (the cutoff-free
physical case), a truncation of the formally exact hierarchy of flow
equations, usually at the level of the four-point vertex, is necessary.

The four-point vertex is directly related to the (momentum resolved) static
spin susceptibility which represents the central outcome of the PFFRG
approach. The onset of magnetic ordering is signaled by a divergence of the
susceptibility along the RG flow (which in a finite system typically reduces
to a finite peak or a kink). Accordingly, non-magnetic (and possibly quantum
spin liquid) phases are characterized by an RG flow that remains smooth down
to the lowest accessible $\Lambda$ scales.

Due to the lack of a small parameter in the purely interacting pseudo-fermion
Hamiltonian, the truncation of the flow equations is an - a priori -
uncontrolled procedure. It can be shown, however, that within the usual
truncation on the level of the four-point vertex, both quantum fluctuations
and classical ordering tendencies are correctly described in leading orders of
$1/N$ and $1/S$, respectively [41, 44]. Here $N$ and $S$ describe the
artificial enlargement of the spin’s symmetry group [$SU(2)\rightarrow SU(N)$]
and the spin length [$1/2\rightarrow S$], respectively. In two very recent
works, certain contributions of the six-point vertex have been taken into
account using a multiloop extension [63, 62] equivalent to a solution of the
parquet self-consistency equations [70, 71, 72]. The quantitative robustness
of the results with respect to increasing loop orders was interpreted as
further evidence for the accuracy of the PFFRG.

##  III $SO(3)$ Majorana Representation

In this section we discuss the $SO(3)$ Majorana representation [64, 65] for
spin-1/2 in detail. For each spin $S^{\alpha}_{i}$ at site $i$, three
different flavors $\alpha\in\\{x,y,z\\}$ of Majorana fermions
$\eta^{\alpha\dagger}_{i}=\eta^{\alpha}_{i}$ are introduced. They fulfill the
anticommutation relations
$\\{\eta^{\alpha}_{i},\eta^{\beta}_{j}\\}=\delta_{ij}\delta^{\alpha\beta}$
which imply $(\eta_{i}^{\alpha})^{2}=1/2$. The formal Hilbert space dimension
per Majorana is $\sqrt{2}$ as appropriate for half a (complex) fermion. The
spin operators
$S^{\alpha}_{i}=-\frac{i}{2}\sum_{\beta\gamma}\varepsilon^{\alpha\beta\gamma}\eta^{\beta}_{i}\eta^{\gamma}_{i}$,
more explicitly written as

$S^{x}_{i}=-i\eta^{y}_{i}\eta^{z}_{i}\text{,}\qquad
S^{y}_{i}=-i\eta^{z}_{i}\eta^{x}_{i}\text{,}\qquad
S^{z}_{i}=-i\eta^{x}_{i}\eta^{y}_{i}\text{,}$ (3)

can be easily checked to fulfill the spin-1/2 algebra

$S_{i}^{\alpha}S_{i}^{\beta}=\frac{1}{4}\delta^{\alpha\beta}+\frac{i}{2}\sum_{\alpha\beta\gamma}\varepsilon^{\alpha\beta\gamma}S_{i}^{\gamma}\text{.}$
(4)

As an example, a Heisenberg coupling term from Hamiltonian (1) is represented
as

$\sum_{\alpha}S^{\alpha}_{i}S^{\alpha}_{j}=-(\eta_{i}^{y}\eta_{i}^{z}\eta_{j}^{y}\eta_{j}^{z}+\eta_{i}^{x}\eta_{i}^{z}\eta_{j}^{x}\eta_{j}^{z}+\eta_{i}^{x}\eta_{i}^{y}\eta_{j}^{x}\eta_{j}^{y}).$
(5)

As usual for auxiliary particle representations, the $SO(3)$ Majorana
representation comes with a gauge freedom. The local $\mathds{Z}_{2}$ gauge
transformation $\eta^{\alpha}_{i}\rightarrow\varepsilon_{i}\eta^{\alpha}_{i}$
with $\varepsilon_{i}=\pm 1$ leaves spin operators invariant since each spin
consists of a product of exactly two Majoranas with equal lattice index. This
gauge freedom is also relevant to understand the structure of the Majorana
Hilbert space. To see this, define the Majorana operator

$\tau_{i}=-2i\eta^{x}_{i}\eta^{y}_{i}\eta^{z}_{i}\text{,}$ (6)

which anticommutes with any $\tau_{j}$ from a different site $j\neq i$ and
fulfills

$\displaystyle\tau_{i}\eta^{\alpha}_{j}=\begin{cases}\eta^{\alpha}_{i}\tau_{i}&\mbox{if
}i=j\\\ -\eta^{\alpha}_{j}\tau_{i}&\mbox{if }i\neq j\\\ \end{cases}\text{.}$
(7)

Consequently, $\tau_{i}$ commutes with all spin operators and thus with any
spin Hamiltonian. To construct a set of mutually commuting operators one needs
to pair $\tau_{i}$ with another conserved Majorana operator.

One choice [73] is to define an additional Majorana $\eta_{i}^{0}$ per site,
so that the parity $p_{i}=2i\tau_{i}\eta_{i}^{0}$ with eigenvalues $\pm 1$ is
a constant of motion. These eigenvalues split the local Majorana Hilbert space
of dimension four into two dynamically decoupled two-dimensional parts each of
which are in one-to-one correspondence to the original local spin Hilbert
space. To invoke $\eta_{i}^{0}$ in the Hamiltonian, parity projection schemes
are required that eventually lead to one of two alternative four-Majorana spin
representations [75]. However, as stated above, we will avoid this additional
complication in the remainder of this work.

An alternative, non-local pairing scheme which does not introduce additional
degrees of freedom requires an even number of sites $N$ [74]. Given an
arbitrary but fixed pairing of sites $(i,j)$, we can define the $N/2$ parities
$p_{(i,j)}=2i\tau_{i}\tau_{j}=\pm 1$. Similar to above, each eigenstate of a
spin Hamiltonian is $2^{N/2}$-fold degenerate, each copy labeled by the above
parities. In other words, the total Majorana Hilbert space dimension of
$2^{3N/2}$ is organized into the usual $2^{N}$ physical spin configurations,
each with an artificial degeneracy of $2^{N/2}$. Choosing a different pairing
of sites corresponds to a unitary rotation of the $2^{N/2}$ basis vectors for
the artificial part of the Hilbert space. Note that since Eq. 3 fully
reproduces the correct spin algebra without the need for an additional
constraint, this Hilbert-space enlargement introduces no unphysical states,
but only exact copies of the physical spin states [75]. This degeneracy is
closely connected to the aforementioned local $\mathds{Z}_{2}$ gauge symmetry:
As the transformation $\tau_{i}\rightarrow-\tau_{i}$ flips the parity
$p_{(i,j)}$, it switches between degenerate states of different parities.

For thermodynamic properties, the above degeneracy leads to the relation
$\mathcal{Z}_{pm}=2^{N/2}\mathcal{Z}$ between the exact partition functions
defined in spin and $SO(3)$ pseudo-Majorana (pm) Hilbert space. Thus, we have
for the physical free energy per site, $f=-T\log\left(\mathcal{Z}\right)/N$,

$\displaystyle f$ $\displaystyle=$ $\displaystyle
f_{pm}+\frac{T}{2}\log\left(2\right)$ (8)

where the first term
$f_{pm}\equiv-\frac{T}{N}\log\left(\mathcal{Z}_{pm}\right)$ will be computed
via PMFRG and the second term accounts for the redundancy inherent in the
$SO(3)$ Majorana representation.

Any expectation values for spin operators (or correlators) $\mathcal{O}_{s}$
are easily computed in the Majorana representation as well. This follows from
the observation that the Majorana version of such an operator,
$\mathcal{O}_{pm}$, is diagonal in the parity sector and the same is true for
any physical density matrix $\rho_{pm}$, like for example the Boltzmann factor
$\rho_{pm}\sim e^{-\beta H_{pm}}$. Then the degeneracy factor $2^{N/2}$ simply
cancels [76] and we have

$\left\langle\mathcal{O}_{s}\right\rangle\equiv\frac{\mathrm{tr}\,\mathcal{O}_{s}\rho_{s}}{\mathrm{tr}\,\rho_{s}}=\frac{\mathrm{tr}\,\mathcal{O}_{pm}\rho_{pm}}{\mathrm{tr}\,\rho_{pm}}\equiv\left\langle\mathcal{O}_{pm}\right\rangle.$
(9)

Finally, we discuss the role of rotations in spin space. In order to employ
the global $SO(3)$ symmetry of the Heisenberg Hamiltonian in Eq. 1 later on,
we demonstrate here that the three Majoranas transform under $SO(3)$ rotations
like the coordinates of a physical vector. Using $\tau_{i}$, the spin
operators can be re-expressed as

$S^{\alpha}_{i}=\tau_{i}\eta^{\alpha}_{i}\text{.}$ (10)

We may now consider the general $SO(3)$ transformation
$\eta^{\alpha}_{i}\rightarrow\sum_{\beta}R_{\alpha\beta}\eta^{\beta}_{i}$ with
$R_{\alpha\beta}\in SO(3)$ being a three dimensional rotation matrix. As
$\tau_{i}$ is invariant under this transformation [75], spin operators must
transform as

$R_{\alpha\beta}S^{\beta}_{i}=\tau_{i}\sum_{\beta}R_{\alpha\beta}\eta^{\beta}_{i}\text{.}$
(11)

It follows that physical $SO(3)$ rotations of a spin $i$ are equivalent to
rotations of the Majorana vector
$\left(\eta^{x}_{i},\eta^{y}_{i},\eta^{z}_{i}\right)$.

##  IV General Majorana FRG Flow Equations

As a basis for our FRG treatment of spin systems in pseudo-Majorana
representation, we first introduce flow equations that are valid for general
interacting Majorana Hamiltonians. To the best of our knowledge, such
equations have not been published in the literature before. We consider

$\displaystyle H=$
$\displaystyle\frac{i}{2}\sum_{\mu_{1,2}}{A_{\mu_{1}\mu_{2}}\eta_{\mu_{1}}\eta_{\mu_{2}}}$
$\displaystyle+$
$\displaystyle\frac{1}{4!}\sum_{\mu_{1,2,3,4}}{V_{\mu_{1}\mu_{2}\mu_{3}\mu_{3}}\eta_{\mu_{1}}\eta_{\mu_{2}}\eta_{\mu_{3}}\eta_{\mu_{4}}}\text{,}$
(12)

where $\\{\mu_{i}\\}$ is an arbitrary set of single-particle indices. As
above, we use the convention
$\left\\{\eta_{\mu_{i}},\eta_{\mu_{j}}\right\\}=\delta_{\mu_{i}\mu_{j}}$.
Majorana exchange statistics require the antisymmetry of $A$ and $V$ under
exchange of any two indices, hermiticity mandates that both couplings must be
real.

Assuming thermal equilibrium, we move on to an imaginary time path integral
formulation [76, 77] defined in terms of Grassmann fields $\eta_{\mu}(\tau)$.
The action reads

$S=\int_{0}^{\beta}{d\tau\left(\sum_{\mu}{\frac{1}{2}\eta_{\mu}(\tau)\partial_{\tau}\eta_{\mu}(\tau)+H\left(\left\\{\eta_{\mu}(\tau)\right\\}\right)}\right)}\text{,}$
(13)

where $\partial_{\tau}$ denotes a derivative with respect to imaginary time
and $\beta=1/T$. We define the Fourier transform
$\eta_{\mu}(\tau)=T\sum_{n}{e^{i\omega_{n}\tau}\eta_{\mu}(i\omega_{n})}$ where
the fermionic Matsubara frequencies are given by $i\omega_{n}=\pi T(2n+1)$,
with $n\in\mathbb{Z}$. In slight abuse of notation, in the following, we will
denote $\omega_{n_{1}}$ by $\omega_{1}$ and equivalently for other
frequencies. The non-interacting part of the action may then be written as

$S_{0}=-\frac{1}{2}\frac{1}{\beta^{2}}\sum_{\omega_{1,2}}{\sum_{\mu_{1,2}}{\eta_{\mu_{1}}(\omega_{1})\left[\bm{G}^{-1}_{0}\right]_{\mu_{1}\omega_{1},\
\mu_{2}\omega_{2}}\eta_{\mu_{2}}(\omega_{2})}}\text{.}$ (14)

with the bare Majorana Green’s function

$\displaystyle\left[\bm{G}^{-1}_{0}\right]_{\mu_{1}\omega_{1},\
\mu_{2}\omega_{2}}$
$\displaystyle=\left(i\omega_{1}\delta_{\mu_{1}\mu_{2}}-iA_{\mu_{1}\mu_{2}}\right)\beta\delta_{\omega_{1},-\omega_{2}}\text{.}$
(15)

This definition is analogous to the complex fermionic bare Green’s function
except for the opposing signs of the two frequencies in the Kronecker delta
related to the absence of an independent Grassmann partner field $\bar{\eta}$
with a relative sign in the Fourier transform.

We are now ready to apply the general FRG scheme from Ref. [17], derived for
an action of a superfield vector $\bm{\Psi}$ containing an arbitrary number of
bosonic or Grassmann fields labeled by the composite index
$l=(\omega_{l},\mu_{l})$,

$\displaystyle S[\bm{\Psi}]$
$\displaystyle=S_{0}[\bm{\Psi}]+S_{\text{int}}[\bm{\Psi}]$
$\displaystyle=-\frac{1}{2}\int_{l}\int_{l^{\prime}}{\Psi_{l}\left[\bm{G}^{-1}_{0}\right]_{l,l^{\prime}}\Psi_{l^{\prime}}}+S_{\text{int}}[\bm{\Psi}]\text{.}$
(16)

where $\int_{l}=\beta^{-1}\sum_{\omega_{l}}\sum_{\mu_{l}}$. A comparison of
Eq. 16 and Eq. 14 yields the direct correspondence
$\bm{\Psi}_{l=(\mu_{l},\omega_{l})}=\eta_{\mu_{l}}(\omega_{l})$. We emphasize
the difference to the superfield vectors of complex fermions or bosons, which
require an additional but independent superfield label, i.e.
$\bm{\Psi}=(\bar{\psi},\psi)$.

The starting point of the FRG scheme is the introduction of a cutoff scale
$\Lambda$ in the bare Green’s function
$\bm{G}_{0}\rightarrow\bm{G}_{0}^{\Lambda}$ such that
$\bm{G}_{0}^{\Lambda=\infty}=0$ and $\bm{G}_{0}^{\Lambda=0}=\bm{G}_{0}$.
Although the flow equations describing the evolution of irreducible vertices
with $\Lambda$ [17] below are general, in the rest of this work, we will
consider a multiplicative Matsubara frequency cutoff
$\Theta^{\Lambda}(\omega_{1})$ to the bare Green’s function

$\displaystyle\left[\bm{G}^{\Lambda}_{0}\right]_{\mu_{1}\omega_{1},\mu_{2}\omega_{2}}$
$\displaystyle=\Theta^{\Lambda}(|\omega_{1}|)\left[\bm{G}_{0}\right]_{\mu_{1}\omega_{1},\mu_{2}\omega_{2}}.$
(17)

At zero temperature, this cutoff is often chosen to be a Heaviside function
$\Theta^{\Lambda}(|\omega|)=\theta(|\omega|-\Lambda)$, at finite temperatures
a smooth cutoff must be chosen instead. While a momentum based cutoff is also
used in some works, we will not consider such schemes here, as our main focus
lies on pseudo-Majoranas without kinetic energy.

As a consequence of the cutoff, the self-energy $\Sigma$ and the four-point
vertex $\Gamma$ acquire implicit dependence on $\Lambda$. These quantities are
defined via the Dyson equation in a superspace spanned by
$(\omega_{i},\mu_{i})$

$\bm{G}=\left[\bm{G}_{0}^{-1}-\bm{\Sigma}\right]^{-1}$ (18)

and the tree-expansion for the connected Green’s functions

$\displaystyle
G^{4,c}_{l_{1},l_{2},l_{3},l_{4}}=-\int_{l_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}}$
$\displaystyle{\bm{G}_{l_{1}l_{1^{\prime}}}\bm{G}_{l_{2}l_{2^{\prime}}}\bm{G}_{l_{3}l_{3^{\prime}}}\bm{G}_{l_{4}l_{4^{\prime}}}}$
$\displaystyle\times\Gamma_{l_{1^{\prime}}l_{2^{\prime}}l_{3^{\prime}}l_{4^{\prime}}}$
(19)

respectively. This $\Lambda$-dependence is given by coupled differential
equations, referred to as flow equations. Physical results can be extracted
from the solution at $\Lambda=0$. Since the action for Majorana systems was
rephrased in superfield notation, we can employ the associated general flow
equations [17] for $\Sigma^{\Lambda}$ and $\Gamma^{\Lambda}$. As appropriate
in thermal equilibrium, and to simplify notation, we employ a modified version
of the Green’s function and vertices with the frequency conserving delta-
function explicitly spelled out,

$\displaystyle\bm{G}_{\mu_{1}\omega_{1},\mu_{2}\omega_{2}}$
$\displaystyle=G_{\mu_{1}\mu_{2}}(\omega_{2})\beta\delta_{\omega_{1},-\omega_{2}}$
(20a) $\displaystyle\bm{\Sigma}_{\mu_{1}\omega_{1},\mu_{2}\omega_{2}}$
$\displaystyle=\Sigma_{\mu_{1}\mu_{2}}(\omega_{1})\beta\delta_{\omega_{1},-\omega_{2}}$
(20b) $\displaystyle\Gamma_{\mu_{1}\omega_{1},\ \mu_{2}\omega_{2},\
\mu_{3}\omega_{3},\ \mu_{4}\omega_{4}}$
$\displaystyle\equiv\Gamma_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$
$\displaystyle\times\beta\delta_{\omega_{1}+\omega_{2}+\omega_{3}+\omega_{4},0}\text{.}$
(20c)

With the above definition, the Dyson equation for fixed frequency indices,
$\bm{G}_{-\omega,\omega}=\left[\left[\bm{G}^{-1}_{0}\right]_{\omega,-\omega}-\bm{\Sigma}_{\omega,-\omega}\right]^{-1}$,
can be written as

$\displaystyle G(\omega)=\left[i\omega-iA-\Sigma(\omega)\right]^{-1}\text{.}$
(21)

The Green’s function and self-energy defined in Eq. (20a) and (20b) fulfill
$G(\omega)=G^{T}(-\omega)$ and $\Sigma(\omega)=\Sigma^{T}(-\omega)$,
respectively.

We also restrict ourselves to the absence of parity symmetry breaking
(expectation values of odd numbers of Majorana operators vanish) and neglect
the contribution from the six-point vertex. The flow equation for the four-
point vertex then separates into three distinct channels, each of which is
characterized by one of the three bosonic transfer frequencies defined as

$\displaystyle s$
$\displaystyle=\omega_{1}+\omega_{2}=-\omega_{3}-\omega_{4},$ $\displaystyle
t$ $\displaystyle=\omega_{1}+\omega_{3}=-\omega_{2}-\omega_{4},$
$\displaystyle u$
$\displaystyle=\omega_{1}+\omega_{4}=-\omega_{2}-\omega_{3}\text{.}$ (22)

The Majorana flow equations for the interaction correction to the free energy,
self energy and the four-point vertex read [17]

$\displaystyle\frac{d}{d\Lambda}F_{\text{int}}^{\Lambda}$
$\displaystyle=\frac{1}{2}\int_{\nu_{1,2,3,4}}T\sum_{\omega^{\prime}}S^{\Lambda}_{\nu_{1}\nu_{2}}(\omega^{\prime})G^{0,\Lambda}_{\nu_{2}\nu_{3}}(-\omega^{\prime})\left[G^{\Lambda}\right]^{-1}_{\nu_{3}\nu_{4}}(-\omega^{\prime})\Sigma^{\Lambda}_{\nu_{4},\nu_{1}}(\omega^{\prime})$
(23a)
$\displaystyle\frac{d}{d\Lambda}\Sigma^{\Lambda}_{\mu_{1},\mu_{2}}(\omega)$
$\displaystyle=-\frac{1}{2}\int_{\nu_{1,2}}T\sum_{\omega^{\prime}}S^{\Lambda}_{\nu_{1}\nu_{2}}(\omega^{\prime})\Gamma^{\Lambda}_{\nu_{1}\nu_{2}\mu_{1}\mu_{2}}(-\omega^{\prime},\omega^{\prime},\omega,-\omega)$
(23b)
$\displaystyle\frac{d}{d\Lambda}\Gamma^{\Lambda}_{\mu_{1},\mu_{2},\mu_{3},\mu_{4}}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$
$\displaystyle=\int_{\nu_{1,2,3,4}}T\sum_{\omega}S^{\Lambda}_{\nu_{1}\nu_{2}}(\omega)$
$\displaystyle\times\bigg{[}$
$\displaystyle\Gamma^{\Lambda}_{\mu_{1}\mu_{2}\nu_{4}\nu_{1}}(\omega_{1},\omega_{2},\omega-s,-\omega)\Gamma^{\Lambda}_{\nu_{2}\nu_{3}\mu_{3}\mu_{4}}(\omega,-\omega+s,\omega_{3},\omega_{4})G^{\Lambda}_{\nu_{3}\nu_{4}}(\omega-s)$
$\displaystyle+$
$\displaystyle\Gamma^{\Lambda}_{\mu_{1}\nu_{1}\mu_{3}\nu_{4}}(\omega_{1},-\omega,\omega_{3},\omega-t)\Gamma^{\Lambda}_{\nu_{2}\mu_{2}\nu_{3}\mu_{4}}(\omega,\omega_{2},-\omega+t,\omega_{4})G^{\Lambda}_{\nu_{3}\nu_{4}}(\omega-t)$
$\displaystyle-$
$\displaystyle\Gamma^{\Lambda}_{\mu_{1}\nu_{4}\nu_{1}\mu_{4}}(\omega_{1},\omega-u,-\omega,\omega_{4})\Gamma^{\Lambda}_{\nu_{3}\mu_{2}\mu_{3}\nu_{2}}(-\omega+u,\omega_{2},\omega_{3},\omega)G^{\Lambda}_{\nu_{3}\nu_{4}}(\omega-u)\bigg{]}\text{.}$
(23c)

As the free energy does not feed back into the other flow equations it is
usually not considered within FRG schemes. In this work, we use its solution
to derive further thermodynamic quantities. In these expressions, we have
introduced the single-scale propagator which is defined as a matrix product of
Green’s functions

$\displaystyle\bm{S}^{\Lambda}$
$\displaystyle=-\bm{G}^{\Lambda}\left[\frac{d}{d\Lambda}\left[\bm{G}^{\Lambda}_{0}\right]^{-1}\right]\bm{G}^{\Lambda}$
$\displaystyle S^{\Lambda}(\omega_{2})$
$\displaystyle=-G^{\Lambda}(\omega_{2})\left[\frac{d}{d\Lambda}\left[G^{\Lambda}_{0}\right]^{-1}(\omega_{2})\right]G^{\Lambda}(\omega_{2})\text{.}$
(24)

In order to solve the flow equations, initial conditions for self-energy and
the four-point vertex are required. As the bare propagator vanishes in this
limit, we immediately see that

$\displaystyle F_{\text{int}}^{\Lambda\rightarrow\infty}$ $\displaystyle=0,$
$\displaystyle\Sigma^{\Lambda\rightarrow\infty}_{\mu_{1},\mu_{2}}(\omega)$
$\displaystyle=0,$
$\displaystyle\Gamma^{\Lambda\rightarrow\infty}_{\mu_{1},\mu_{2},\mu_{3},\mu_{4}}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$
$\displaystyle=V_{\mu_{1},\mu_{2},\mu_{3},\mu_{4}}\text{.}$ (25)

##  V Symmetry-Based Vertex Parametrization

We now specialize the general Majorana FRG of this section to treat the
interacting system of pseudo-Majoranas ensuing from the application of the
representation (3) to the Heisenberg spin-1/2 Hamiltonian (1),

$H=-\sum_{(i,j)}J_{ij}{\left(\eta_{i}^{y}\eta_{j}^{y}\eta_{i}^{z}\eta_{j}^{z}+\eta_{i}^{z}\eta_{j}^{z}\eta_{i}^{x}\eta_{j}^{x}+\eta_{i}^{x}\eta_{j}^{x}\eta_{i}^{y}\eta_{j}^{y}\right)}.$
(26)

As a first step, we proceed with a detailed discussion of the parametrization
of vertices and propagators using the symmetries of our model. Following the
approach of Ref. [54], we will first derive symmetry relations for the Green’s
functions defined as

$\displaystyle G(1,2)$
$\displaystyle=\int_{0}^{\beta}{d\tau_{1}d\tau_{2}e^{i\omega_{1}\tau_{1}}e^{i\omega_{2}\tau_{2}}\expectationvalue{\eta_{\mu_{1}}(\tau_{1})\eta_{\mu_{2}}(\tau_{2})}}$
$\displaystyle=G_{\mu_{1},\mu_{2}}(\omega_{2})\beta\delta_{\omega_{1},-\omega_{2}}$
(27) $\displaystyle G^{4}(1,2,3,4)$
$\displaystyle=\int_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}d\tau_{4}e^{i(\omega_{1}\tau_{1}+\omega_{2}\tau_{2}+\omega_{3}\tau_{3}+\omega_{4}\tau_{4})}$
$\displaystyle\times\expectationvalue{\eta_{\mu_{1}}(\tau_{1})\eta_{\mu_{2}}(\tau_{2})\eta_{\mu_{3}}(\tau_{3})\eta_{\mu_{4}}(\tau_{4})}$
(28)
$\displaystyle=G^{4}_{\mu_{1},\mu_{2},\mu_{3},\mu_{4}}(s,t,u)\beta\delta_{\omega_{1}+\omega_{2}+\omega_{3}+\omega_{4},0}\text{.}$
(29)

where the labels $(1,2,3,4)$ contain all arguments that are not explicitly
specified, i.e $1=(\mu_{1},\omega_{1})$ in this case. Matsubara frequency
conservation follows from the fact that thermal expectation values only depend
on imaginary time differences. The time-ordering operator is suppressed since
it is included in the path integral formalism by default. The properties
derived in the following will then carry over to $\Sigma$ and $\Gamma$ due to
their relations via Eqs. 18 and 19.

###  V.1 Hermiticity

The Hamiltonian is a hermitian operator, satisfying $H=H^{\dagger}$. Due to
$\expectationvalue{\mathcal{O}}^{*}=\expectationvalue{\mathcal{O^{\dagger}}}$
and $\eta(\tau)^{\dagger}=e^{-H\tau}\eta e^{H\tau}=\eta(-\tau)$ in the
Heisenberg picture, one can find the complex conjugate of the two-point
Green’s functions as

$\displaystyle G(1,2)^{*}$
$\displaystyle=\int{d\tau_{1}d\tau_{2}e^{-i\omega_{1}\tau_{1}-i\omega_{2}\tau_{2}}\expectationvalue{\eta_{\mu_{2}}(-\tau_{2})\eta_{\mu_{1}}(-\tau_{1})}}$
$\displaystyle=-G(1,2)\text{.}$ (30)

As a consequence, the two-point Green’s function in Matsubara frequency space
is purely imaginary and from an analogous argument, the four-point Green’s
function must be real

$\displaystyle G(1,2)\in i\mathds{R},$ $\displaystyle
G^{4}(1,2,3,4)\in\mathds{R}\text{.}$ (31)

###  V.2 Time reversal symmetry

Time reversal $T$ is an anti-unitary operation
$(\left\langle\psi|\psi^{\prime}\right\rangle^{*}=\left\langle
T\psi|T\psi^{\prime}\right\rangle)$ which in the present case can be defined
by performing a complex conjugation while leaving Majorana operators invariant
[78]:

$TiT^{-1}=-i\text{,}\>\>\quad T\eta_{\mu}T^{-1}=\eta_{\mu}\text{.}$ (32)

This flips the sign of the spin operators (3) as required. Time reversal
symmetry is violated by an external magnetic field or, more generally, any
Majorana bilinear in the Hamiltonian. For a $T$-symmetric Hamiltonian
$THT^{-1}=H$, thermal expectation values obey
$\left\langle\mathcal{O}\right\rangle=\left\langle
T\mathcal{O}T^{-1}\right\rangle^{*}$. From this, we have
$\expectationvalue{\eta_{\mu_{1}}(\tau_{1})\eta_{\mu_{2}}(\tau_{2})}=\expectationvalue{\eta_{\mu_{1}}(\tau_{1})\eta_{\mu_{2}}(\tau_{2})}^{*}$
and with Eq. (27), it follows that

$G_{\mu_{1}\mu_{2}}(\omega_{1},\omega_{2})=G_{\mu_{1}\mu_{2}}(-\omega_{1},-\omega_{2})^{*}\text{.}$
(33)

Similarly, the four-point correlator has the property

$G^{4}(1,2,3,4)=G^{4}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(-\omega_{1},-\omega_{2},-\omega_{3},-\omega_{4})^{*}\text{.}$
(34)

###  V.3 Local $\mathds{Z}_{2}$ gauge redundancy

Since our considerations from here on require the explicit specification of
site indices, we will now separate the previously used superlabel $\mu$ into a
site-index and a Majorana flavor $\mu\rightarrow(i,\alpha)$. In the $SO(3)$
Majorana representation spins are invariant under the gauge transformation
$\eta^{\alpha}_{i}\rightarrow\varepsilon_{i}\eta^{\alpha}_{i}$ for all
$\alpha=x,y,z$ with $\varepsilon_{i}=\pm 1$ for an arbitrary lattice site $i$.
Since expectation values must be invariant under gauge transformations as
well, we may write

$\expectationvalue{\eta^{\alpha_{1}}_{i_{1}}(\tau_{1})\eta^{\alpha_{2}}_{i_{2}}(\tau_{2})}=\varepsilon_{i_{1}}\varepsilon_{i_{2}}\expectationvalue{\eta^{\alpha_{1}}_{i_{1}}(\tau_{1})\eta^{\alpha_{2}}_{i_{2}}(\tau_{2})}\text{,}$
(35)

where $\varepsilon_{i_{1}}\varepsilon_{i_{2}}=-1$ may always be chosen for two
different sites. As a consequence, non-zero propagators must contain an even
number of Majorana operators from each site, so that

$G_{i_{1}i_{2}}(1,2)\equiv\delta_{i_{1}i_{2}}G_{i_{1}}(1,2)\text{.}$ (36)

Likewise, the four-point correlator can only depend on up to two distinct
sites only, so we choose

$G^{4}_{i,i,j,j}(1,2,3,4)\equiv G^{4}_{ij}(1,2;3,4)\text{.}$ (37)

Correlators of the form $ijij$ and $ijji$ need to be brought to the standard
form Eq. 37 using fermionic anticommutation rules, which restricts the number
of allowed permutations in $G^{4}_{ij}(1,2;3,4)$ to exchanges of the first and
last two indices only. As a consequence of the (bi-)local nature of
propagators (four-point vertices), the site summations in the flow equations
can be simplified. The special case $i=j$ for the four-point vertex needs to
be considered separately. The corresponding flow equations can then be
expressed diagrammatically as shown in Fig. 1. The bubble-diagram
corresponding to the $s$-channel of the non-local vertex $\Gamma_{ij}$ shown
in Fig. 1 d) is of particular interest. As in the PFFRG this diagram includes
the random-phase approximation which is responsible for the emergence of long-
range magnetic order [41].

Figure 1: $\mathds{Z}_{2}$-invariant Majorana FRG flow equations for the
interaction correction to the free energy (a), the self-energy (b), and the
local (c) and nonlocal (d) four-point vertices. The order of labels
$1=(\alpha_{1},\omega_{1})$ always corresponds to that on the left hand side
of the vertex flow equations such that the site index is conserved along solid
lines. In these equations, internal lines correspond to fully dressed Green’s
functions $G_{i}(1,2)$, while the single scale propagator $S_{i}(1,2)$ is
represented by a slashed line. Similarly, the crossed line in a) corresponds
to the local propagator $\left[SG^{0}G^{-1}\right]_{i}(1,2)$.

###  V.4 Lattice symmetries

For simplicity, the systems that are considered in the following consist of
equivalent sites. Correlators can then always be computed with one arbitrary
reference site fixed. Combining this with local $\mathds{Z}_{2}$ gauge
redundancy eliminates all site indices of the two-point correlator. Similarly,
four-point correlators depend only on the distance vector between the two
sites. Although this means that the order of site indices in $\Gamma_{ij}$ is
irrelevant for systems with equivalent sites, we will not make use of this
property. As a result, the pseudo-Majorana flow equations presented here are
easily generalized towards non-Bravais lattices by adding an additional
sublattice-index. Most lattice systems further exhibit point-group symmetries,
such as the $C_{4}$ rotation symmetry and mirror planes of the square lattice,
which can straightforwardly be used to reduce the numerical effort and are not
further discussed in the following due to their lattice-specific nature.

###  V.5 Global $SO(3)$ rotation symmetry

The global $SO(3)$ spin-rotation symmetry of the Heisenberg model can easily
be translated to vertex functions. As discussed in Sec. III, global spin
rotations specified by a $3\times 3$ rotation matrix
$R_{\alpha\mu}(\bm{\phi})$ act on the Majorana fermions as

$\eta^{\alpha}_{i}\rightarrow\sum_{\beta}R_{\alpha\beta}(\bm{\phi})\eta^{\beta}_{i}\quad\forall
i\text{.}$ (38)

The Heisenberg Hamiltonian is invariant under spin rotations due to the
isotropic nature of its couplings.

We will now apply this symmetry to restrict the types of vertices and find
relations between vertices with different flavor indices. Of particular
interest are the specific rotations along the $x$, $y$ and $z$-axes as
displayed in Table 1. The combination $R_{x}(\pi/2)\circ
R_{z}(\pi/2)\equiv\mathcal{P}$ realizes an anti-cyclic permutation of the
flavors.

Angle | x | y | z
---|---|---|---
$\pi/2$ | $\eta^{y}\rightarrow-\eta^{z}$ | $\eta^{x}\rightarrow\eta^{z}$ | $\eta^{x}\rightarrow-\eta^{y}$
| $\eta^{z}\rightarrow\eta^{y}$ | $\eta^{z}\rightarrow-\eta^{x}$ | $\eta^{y}\rightarrow\eta^{x}$
$\pi$ | $\eta^{y}\rightarrow-\eta^{y}$ | $\eta^{x}\rightarrow-\eta^{x}$ | $\eta^{x}\rightarrow-\eta^{x}$
| $\eta^{z}\rightarrow-\eta^{z}$ | $\eta^{z}\rightarrow-\eta^{z}$ | $\eta^{y}\rightarrow-\eta^{y}$
Table 1: Symmetry transformations corresponding to specific spin rotations
along the $x,y$ and $z$ axes.

We apply these symmetries to correlators, using the convention
$\gamma\neq\alpha\neq\beta\neq\gamma$ to refer to fixed, pairwise different
flavors. In this way, we find that the two-point Green’s function does not
depend on any flavor labels.

$\displaystyle\expectationvalue{\eta^{\alpha}_{1}\eta^{\beta}_{2}}$
$\displaystyle\stackrel{{\scriptstyle
R_{\alpha}(\pi)}}{{=}}-\expectationvalue{\eta^{\alpha}_{1}\eta^{\beta}_{2}}=0,$
$\displaystyle\Rightarrow G_{\alpha_{1},\alpha_{2}}(1,2)$
$\displaystyle=G_{\alpha_{1}}(1,2)\delta_{\alpha_{1},\alpha_{2}}\stackrel{{\scriptstyle\mathcal{P}}}{{=}}G(1,2)\delta_{\alpha_{1},\alpha_{2}}\text{.}$
(39)

Because the four-point correlator has four flavor indices, at least two of
them must be equal. An argument analogous to above shows that only vertices
with an even number of flavors can be nonzero. Furthermore, rotations by
$\pi/2$ transform different flavor combinations into each other, for instance

$\expectationvalue{\eta^{\alpha}_{1}\eta^{\alpha}_{2}\eta^{\beta}_{3}\eta^{\beta}_{4}}\stackrel{{\scriptstyle
R_{\alpha}(\pi/2)}}{{=}}\expectationvalue{\eta^{\alpha}_{1}\eta^{\alpha}_{2}\eta^{\gamma}_{3}\eta^{\gamma}_{4}}\text{.}$
(40)

These arguments identify four independent flavor configurations for the four-
point correlator, $G^{4}_{xxxx}(1,2,3,4)$, $G^{4}_{xxyy}(1,2,3,4)$,
$G^{4}_{xyxy}(1,2,3,4)$ and $G^{4}_{xyyx}(1,2,3,4)$, all other types are
either zero or related by Eq. 40.

After these simplifications, we consider a general rotation to find a relation
between those four different correlators. Since they are now parametrized in
terms of $x$ and $y$, we only need to consider rotations along the $z$-axis.
The $\eta^{x}$ Majoranas then transform as
$\eta^{x}_{i}\rightarrow\cos\theta\eta^{x}_{i}-\sin\theta\eta^{y}_{i}$ so that

$G^{4}_{xxxx}\stackrel{{\scriptstyle
R_{z}(\theta)}}{{=}}\expectationvalue{\left(\cos\theta\eta^{x}_{1}-\sin\theta\eta^{y}_{1}\right)\dots\left(\cos\theta\eta^{x}_{4}-\sin\theta\eta^{y}_{4}\right)}\text{.}$
(41)

Expanding the product and using the above symmetries, we obtain a relation
independent of $\theta$

$G^{4}_{xxxx}=G^{4}_{xxyy}+G^{4}_{xyxy}+G^{4}_{xyyx}\text{,}$ (42)

where the argument $(1,2,3,4)$ has been suppressed. Since we considered an
arbitrary rotation, our last consideration further serves as a proof that no
other symmetries than the ones already shown may be found from $SO(3)$
rotations. Indeed, one arrives at the same identity regardless of which type
of correlator one transforms (i.e. transforming $G^{4}_{xyxy}$ yields the same
result). Rotations along the $x$ or $y$ direction also generate no further
information as a result of the permutation symmetry $\mathcal{P}$ and
rotations around an arbitrary axis may always be decomposed as a product of
$x,y$ and $z$ rotations.

##  VI Pseudo-Majorana fRG Flow Equations

The symmetries of the last section imply the following parametrization of the
pseudo-Majorana propagator,

$G(1,2)=\mathcal{G}(-\omega_{1})\delta_{i_{1},i_{2}}\delta_{\alpha_{1},\alpha_{2}}\delta_{\omega_{1},-\omega_{2}}\beta\text{,}$
(43)

where the imaginary and antisymmetric self-energy, abbreviated as
$\Sigma(\omega)=-i\gamma(\omega)$, enters via the Dyson equation (18),

$\mathcal{G}(\omega)=\frac{1}{i\omega+\ i\gamma(\omega)}\equiv-ig(\omega).$
(44)

In analogy to the real functions $\gamma(\omega)$ and $g(\omega)$ we also
replace the imaginary single scale propagator via
$S^{\Lambda}(\omega)=-i\dot{g}^{\Lambda}(\omega)$. Due to the diagonal
structure of the propagators, the symmetries for the four-point Green’s
functions then carry over to vertex functions [cf. Eq. 19] whose frequency
dependence is parametrized by the three bosonic frequencies introduced in Eq.
22. The three independent four-point vertices are

$\displaystyle\Gamma_{a\ ij}(s,t,u)$ $\displaystyle\equiv\Gamma_{xi,\ xi,\
xj,\ xj}(s,t,u),$ $\displaystyle\Gamma_{b\ ij}(s,t,u)$
$\displaystyle\equiv\Gamma_{xi,\ xi,\ yj,\ yj}(s,t,u),$
$\displaystyle\Gamma_{c\ ij}(s,t,u)$ $\displaystyle\equiv\Gamma_{xi,\ yi,\
xj,\ yj}(s,t,u)\text{.}$ (45)

In the special case $i=j$, there are only two independent vertices since

$\Gamma_{c\ ii}(s,t,u)=-\Gamma_{b\ ii}(t,s,u)\text{.}$ (46)

Vertices with negative bosonic frequencies are symmetry related to positive
frequencies by time-reversal and a symmetry $t\leftrightarrow u$ further
allows to reduce the numerical effort. Details are given in Table 2.

Operation | Symmetry for $\Gamma_{\mu\ ij}(s,t,u)$ | valid $\mu$
---|---|---
$1\leftrightarrow 2$ | $t\leftrightarrow u\text{\ and\ }\Gamma_{\mu}\leftrightarrow-\Gamma_{\mu}$ | $a,b$
$T\circ(1,3)\leftrightarrow(2,4)$ | $s\leftrightarrow-s$ | $a,b,c$
$T\circ(1,2)\leftrightarrow(3,4)$ | $t\leftrightarrow-t\text{\ and\ }i\leftrightarrow j$ | $a,b,c$
$T\circ(1,2)\leftrightarrow(4,3)$ | $u\leftrightarrow-u\text{\ and\ }i\leftrightarrow j$ | $a,b,c$
Table 2: Transformations of the frequency arguments under time reversal $T$
and specific permutations of indices in$\Gamma_{ij}(1,2;3,4)$. The latter
three rows apply to all three types of vertices and allow for a
parametrization using positive frequencies only. Note that the final two
permutations also exchange the order of $i$ and $j$ which is of importance for
non-Bravais lattices. The remaining $t\leftrightarrow u$ symmetry for
$\Gamma_{c}$ can be established by the exchange $1\leftrightarrow 2$, which
changes the vertex to the form $\Gamma_{xyyx}$. Using Eq. 42 to express
$\Gamma_{xyyx}(s,t,u)=-\Gamma_{c}(s,u,t)$ in terms of the other vertices used
in the parametrization, we obtain $\Gamma_{c\ ij}(s,u,t)=(-\Gamma_{a\
ij}+\Gamma_{b\ ij}+\Gamma_{c\ ij})(s,t,u)$.

In the above parametrization, the flow equations for the interaction
correction to the free energy per spin and the self-energy may be simplified.
Specifying the external flavor and site indices on the left hand side of the
flow equations, we directly perform flavor sums to obtain

$\frac{d}{d\Lambda}f_{\text{int}}^{\Lambda}=-\frac{3T}{2}\sum_{\omega}\dot{g}^{\Lambda}(\omega)\frac{g^{0,\Lambda}(\omega)}{g^{\Lambda}(\omega)}\gamma^{\Lambda}(\omega),$
(47)
$\displaystyle\frac{d}{d\Lambda}\gamma^{\Lambda}(\omega_{1})=\frac{T}{2}\sum_{\omega}\sum_{j}\dot{g}^{\Lambda}(\omega)\bigg{\\{}$
$\displaystyle\Gamma^{\Lambda}_{a\ ij}(0,\omega_{1}+\omega,\omega_{1}-\omega)$
$\displaystyle+2$ $\displaystyle\Gamma^{\Lambda}_{b\
ij}(0,\omega_{1}+\omega,\omega_{1}-\omega)\bigg{\\}}\text{.}$ (48)

Similarly, we may now express the flow equations for four-point vertices in
the same way. For conciseness of notation, both the initial fermionic
frequencies as well as the exchange frequencies $s,t$ and $u$ will be used on
the right hand side which are defined by Eq. 22, or inversely,

$\displaystyle\omega_{1}$
$\displaystyle=\frac{s+t+u}{2}\text{,}\quad\omega_{2}=\frac{s-t-u}{2}$
$\displaystyle\omega_{3}$
$\displaystyle=\frac{-s+t-u}{2}\text{,}\quad\omega_{4}=\frac{-s-t+u}{2}\text{.}$
(49) $\displaystyle\frac{d}{d\Lambda}\Gamma^{\Lambda}_{a\ ij}(s,t,u)$
$\displaystyle=X^{\Lambda}_{a\ ij}(s,t,u)-\tilde{X}^{\Lambda}_{a\
ij}(t,s,u)+\tilde{X}^{\Lambda}_{a\ ij}(u,s,t)$ (50a)
$\displaystyle\frac{d}{d\Lambda}\Gamma^{\Lambda}_{b\ ij}(s,t,u)$
$\displaystyle=X^{\Lambda}_{b\ ij}(s,t,u)-\tilde{X}^{\Lambda}_{c\
ij}(t,s,u)+\tilde{X}^{\Lambda}_{c\ ij}(u,s,t)$ (50b)
$\displaystyle\frac{d}{d\Lambda}\Gamma^{\Lambda}_{c\ i,j\neq i}(s,t,u)$
$\displaystyle=X^{\Lambda}_{c\ ij}(s,t,u)-\tilde{X}^{\Lambda}_{b\
ij}(t,s,u)+\tilde{X}^{\Lambda}_{d\ ij}(u,s,t)$ (50c) $\displaystyle
X^{\Lambda}_{a\ ij}(s,t,u)$
$\displaystyle=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\sum_{k}\left[\Gamma^{\Lambda}_{a\
ki}\left(s,\omega+\omega_{1},\omega+\omega_{2}\right)\Gamma^{\Lambda}_{a\
kj}\left(s,\omega-\omega_{3},\omega-\omega_{4}\right)+2(a\rightarrow
b)\right]$ (51a) $\displaystyle X^{\Lambda}_{b\ ij}(s,t,u)$
$\displaystyle=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\sum_{k}\left[\Gamma^{\Lambda}_{a\
ki}\left(s,\omega+\omega_{1},\omega+\omega_{2}\right)\Gamma^{\Lambda}_{b\
kj}\left(s,\omega-\omega_{3},\omega-\omega_{4}\right)+(a\rightarrow
b)+(a\leftrightarrow b)\right]$ (51b) $\displaystyle X^{\Lambda}_{c\
ij}(s,t,u)$
$\displaystyle=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\sum_{k}\left[\Gamma^{\Lambda}_{c\
ki}\left(s,\omega+\omega_{1},\omega+\omega_{2}\right)\Gamma^{\Lambda}_{c\
kj}\left(s,\omega-\omega_{3},\omega-\omega_{4}\right)+(\omega_{1}\leftrightarrow\omega_{2},\omega_{3}\leftrightarrow\omega_{4})\right]$
(51c) $\displaystyle\tilde{X}^{\Lambda}_{a\ i,j\neq
i}(s,t,u)=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\big{\\{}\big{[}$
$\displaystyle\Gamma^{\Lambda}_{a\
ij}\left(\omega+\omega_{2},s,\omega+\omega_{1}\right)\Gamma^{\Lambda}_{a\
ij}\left(\omega-\omega_{4},s,\omega-\omega_{3}\right)$ $\displaystyle+$
$\displaystyle(\omega_{1}\leftrightarrow\omega_{2},\omega_{3}\leftrightarrow\omega_{4},i\leftrightarrow
j)\big{]}+2(a\rightarrow c)\big{\\}}$ (52a)
$\displaystyle\tilde{X}^{\Lambda}_{b\ i,j\neq
i}(s,t,u)=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\big{\\{}\big{[}$
$\displaystyle\Gamma^{\Lambda}_{a\
ij}\left(\omega+\omega_{2},s,\omega+\omega_{1}\right)\Gamma^{\Lambda}_{c\
ij}\left(\omega-\omega_{4},s,\omega-\omega_{3}\right)$ $\displaystyle+$
$\displaystyle(\omega_{1}\leftrightarrow\omega_{2},\omega_{3}\leftrightarrow\omega_{4},i\leftrightarrow
j)\big{]}+(a\rightarrow c)+(a\leftrightarrow c)\big{\\}}$ (52b)
$\displaystyle\tilde{X}^{\Lambda}_{c\ i,j\neq
i}(s,t,u)=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\big{\\{}\big{[}$
$\displaystyle\Gamma^{\Lambda}_{b\
ij}\left(\omega+\omega_{2},\omega+\omega_{1},s\right)\Gamma^{\Lambda}_{b\
ij}\left(\omega-\omega_{4},\omega-\omega_{3},s\right)$ $\displaystyle+$
$\displaystyle(\omega_{1}\leftrightarrow\omega_{2},\omega_{3}\leftrightarrow\omega_{4},i\leftrightarrow
j)\big{]}+(b\rightarrow c)\big{\\}}$ (52c)
$\displaystyle\tilde{X}^{\Lambda}_{d\ i,j\neq
i}(s,t,u)=T\sum_{\omega}\dot{g}^{\Lambda}(\omega)g^{\Lambda}(\omega+s)\big{\\{}\big{[}$
$\displaystyle\Gamma^{\Lambda}_{b\
ij}\left(\omega+\omega_{2},\omega+\omega_{1},s\right)\Gamma^{\Lambda}_{c\
ij}\left(\omega-\omega_{4},\omega-\omega_{3},s\right)$ $\displaystyle+$
$\displaystyle(\omega_{1}\leftrightarrow\omega_{2},\omega_{3}\leftrightarrow\omega_{4},i\leftrightarrow
j)\big{]}+(b\leftrightarrow c)\big{\\}}$ (52d)

To reduce the length of expressions, we have defined the single-channel
contributions $X^{\Lambda}_{a,b,c\ ij}$ and $\tilde{X}^{\Lambda}_{a,b,c,d\
ij}$ in Eqs. 51 and 52 [46]. The flow equations of local vertices are obtained
noting that $\tilde{X}^{\Lambda}_{a,b,c\ ii}(s,t,u)\equiv X^{\Lambda}_{a,b,c\
ii}(s,t,u)$. We further stress that no flow equation for $\Gamma_{c\ ii}$ is
required in Eq. 50, as this vertex is equivalent to $\Gamma_{b\ ii}$ by virtue
of Eq. 46.

In the PFFRG, the Katanin truncation scheme [79] was instrumental in providing
sufficient feedback of the self-energy flow into the vertex flow equations
[22]. It amounts to promoting the single-scale propagator in the flow
equations of four-point vertices to a full derivative of the Green’s function

$\displaystyle S^{\Lambda}(\omega)$
$\displaystyle\rightarrow\frac{d}{d\Lambda}\mathcal{G}^{\Lambda}(\omega)\equiv
S_{\text{conv.}}(\omega)+S_{\text{Kat}}(\omega)$
$\displaystyle=-\mathcal{G}(\omega)^{2}\frac{d}{d\Lambda}\left[\mathcal{G}^{0\Lambda}(\omega)\right]^{-1}+\mathcal{G}(\omega)^{2}\frac{d}{d\Lambda}\Sigma^{\Lambda}(\omega)\text{.}$
(53)

At zero temperature, frequencies become continuous and
$T\sum_{\omega}\rightarrow(2\pi)^{-1}\int d\omega$. Using the sharp frequency
cutoff
$\mathcal{G}^{0\Lambda}(\omega)=\mathcal{G}^{0}(\omega)\theta(|\omega|-\Lambda)$,
we thus obtain in the usual way using Morris’s Lemma [80]

$\dot{g}^{\Lambda}_{T=0}(\omega)=-\frac{\delta(|\omega|-\Lambda)}{\omega+\gamma^{\Lambda}(\omega)}+\dot{g}^{\Lambda}_{\text{Kat}}(\omega)\text{.}$
(54)

At finite temperatures, a sharp cutoff of frequencies is no longer possible
due to ambiguities that arise if $|\omega|-\Lambda$ lies between two discrete
Matsubara frequencies. Noting that there is still freedom in the choice of a
smooth cutoff [57, 81], here we choose a Lorentzian cutoff function

$\Theta^{\Lambda}(\omega_{n})=\frac{\omega_{n}^{2}}{\omega_{n}^{2}+\Lambda^{2}}\text{.}$
(55)

Using Eqs. 18, 17 and 53 the expressions for the Green’s function and the
single-scale propagator become

$\displaystyle g^{\Lambda}(i\omega_{n})$
$\displaystyle=\frac{\omega_{n}}{\omega_{n}^{2}+\omega_{n}\gamma(\omega_{n})+\Lambda^{2}}$
$\displaystyle\dot{g}^{\Lambda}(i\omega_{n})$
$\displaystyle=-g^{2}(i\omega_{n})\left(\frac{2\Lambda}{\omega_{n}}+\frac{d\gamma^{\Lambda}(i\omega_{n})}{d\Lambda}\right)\text{.}$
(56)

Finally, we need to specify the initial conditions for the newly defined
vertices. After re-expressing the Heisenberg Hamiltonian (1) by insertion of
Eq. 3 for the spin operators, a comparison of coefficients yields

$\displaystyle f_{\text{int}}^{\Lambda\rightarrow\infty}$ $\displaystyle=0,$
$\displaystyle\Sigma^{\Lambda\rightarrow\infty}$ $\displaystyle=0,$
$\displaystyle\Gamma_{a\ ij}^{\Lambda\rightarrow\infty}$
$\displaystyle=\Gamma_{b\ ij}^{\Lambda\rightarrow\infty}=0,$
$\displaystyle\Gamma_{c\ ij}^{\Lambda\rightarrow\infty}$
$\displaystyle=-J_{ij}\text{.}$ (57)

To summarize, in our PMFRG scheme the flow equations for the free energy (47),
self-energy (48) and the vertex functions (50), are solved numerically
starting from large but finite $\Lambda\gg J$ down to $\Lambda\simeq 0$,
approximating the initial conditions with the $\Lambda\rightarrow\infty$
values presented above. The flow of the free energy correction is integrated
along the way but does not feed back into the other flow equations. Further
details on the numerical implementation of the PMFRG are given in Appendix A.
The next section describes how to extract observables along the flow and, most
importantly, at the physical endpoint $\Lambda=0$.

##  VII Observables

In this section, we discuss the observables for Heisenberg spin-1/2 systems
that will be studied in the following sections. These are the free energy,
internal energy, heat capacity and static susceptibility. We explain how these
observables are calculated from the eigenstates and -energies of the spin
Hamiltonian (1), its exact representation with $SO(3)$ Majorana fermions and
from the (approximate) solution of the PMFRG flow equations.

From the partition function of a $N$-spin system with eigenenergies $E_{n}$,
$\mathcal{Z}=\sum_{n}e^{-\beta E_{n}}$, the free energy per spin is given by

$F/N=f=-\frac{T}{N}\log\left(\mathcal{Z}\right)=-\frac{T}{N}\log\sum_{n}e^{-\beta
E_{n}}.$ (58)

The energy per spin is

$E/N=-\frac{\partial\log\left(\mathcal{Z}\right)}{N\partial\beta}=\frac{\partial(f\beta)}{\partial\beta}=\frac{1}{N\mathcal{Z}}\sum_{n}E_{n}e^{-\beta
E_{n}},$ (59)

which as a function of $T$ also determines the heat capacity

$C/N=\frac{\partial}{\partial
T}E/N=\frac{1}{NT^{2}}\left(\frac{1}{Z}\sum_{n}E_{n}^{2}e^{-\beta
E_{n}}-E^{2}\right).$ (60)

For small systems amenable to exact diagonalization, the rightmost expressions
are most convenient. From the solution of the PMFRG flow equation (47) for the
interaction correction to the pseudo-Majorana free energy per site, we find
$f_{pm}=f_{pm,0}+f_{\text{int}}^{\Lambda=0}$. The non-interacting free energy
for three pseudo-Majoranas per site is $f_{pm,0}=-3T\,\mathrm{log}(2)/2$.
Using the relation between $f_{pm}$ and $f$, Eq. (8), we finally obtain

$f=-T\,\mathrm{log}(2)+f_{\text{int}}^{\Lambda=0}.$ (61)

The static spin-spin correlator can be computed from

$\chi_{ij}=\int_{0}^{\beta}d\tau\expectationvalue{S^{z}_{i}(\tau)S^{z}_{j}(0)}.$
(62)

Note that $\chi_{ij}$ can also be interpreted as a static (zero-field)
susceptibility as it measures the response of a spin at site $i$ when a
magnetic field is exerted at site $j$. We represent the spin operators by
Majorana fermions and obtain from the vertices of the PMFRG at cutoff scale
$\Lambda$,

$\displaystyle\chi^{\Lambda}_{ij}=$
$\displaystyle+T^{2}\sum_{\omega_{1}\omega_{2}}g^{\Lambda}(\omega_{1})^{2}g^{\Lambda}(\omega_{2})^{2}\Gamma^{\Lambda}_{c\
ij}(0,\omega_{1}+\omega_{2},\omega_{1}-\omega_{2})$
$\displaystyle+T\sum_{\omega_{1}}g^{\Lambda}(\omega_{1})^{2}\delta_{ij}\text{.}$
(63)

Of particular interest for the two-dimensional systems below is the uniform
susceptibility $\chi=\sum_{i,j}\chi_{ij}$.

##  VIII Application: Small Spin Clusters

Figure 2: Zero temperature PMFRG flow of the static local and nonlocal
susceptibilities $\chi_{ij}$ for the antiferromagnetic Heisenberg dimer. The
grey dotted line represents the exact physical ($\Lambda=0$) result.

###  VIII.1 Spin dimer and the fermion parity issue

Small spin clusters constitute an ideal testbed for probing the accuracy of
our approaches as they already represent non-trivial problems within the PMFRG
(and PFFRG) but are still exactly solvable. We first investigate the simple
case of two spins, $i=0,1$ coupled with an antiferromagnetic Heisenberg
interaction $J=1$. Due to the small Hilbert space, this dimer model
$H_{N=2}=\sum_{\alpha}S_{0}^{\alpha}S_{1}^{\alpha}$ is analytically solvable.
While the free energy Eq. 58 is straightforwardly found, some care is required
for the calculation of the susceptibility from the Lehmann representation
where the term contributing in the case $i\nu+E_{n}-E_{m}=0$ is often
neglected in textbook derivations. We obtain

$\displaystyle\chi_{00}$
$\displaystyle=\frac{e^{\beta}-1+\beta}{2(e^{\beta}+3)},$
$\displaystyle\chi_{01}$
$\displaystyle=-\frac{e^{\beta}-1-\beta}{2(e^{\beta}+3)}\text{.}$ (64)

Our PMFRG results for the static susceptibility in the case $T=0$ are shown in
Fig. 2 as a function of the cutoff. We find that $\chi_{ij}^{\Lambda}$ flows
smoothly without any feature, surpasses the exact results $\chi_{ij}=\pm 0.5$
and diverges at $\Lambda=0$. This unphysical divergence is not restricted to
the Heisenberg dimer but appears in all other models considered here. However,
the dimer allows for the most simple discussion of the origin of this
divergence, which equally plagues the flow of the nonlocal vertices of type
$\Gamma_{a,01}=\Gamma_{x0,x0,x1,x1}$ and $\Gamma_{c,01}=\Gamma_{x0,y0,x1,y1}$.

To explain the origin of this divergence, consider the Heisenberg dimer which
can be exactly solved in the $SO(3)$ Majorana representation,

$H_{N=2}=-\frac{1}{4}p^{x}p^{y}p^{z}\left(p^{x}+p^{y}+p^{z}\right)\text{.}$
(65)

Here, $p^{\alpha}\equiv 2i\eta^{\alpha}_{0}\eta^{\alpha}_{1}$ are the three
flavor parities related to the non-local parity introduced in Sec. III via
$p_{(0,1)}=2i\tau_{0}\tau_{1}=-p^{x}p^{y}p^{z}$. While $p_{(i,j)}=\pm 1$ is
always conserved for generic spin systems, $p^{\alpha}=\pm 1$ are additional
constants of motion only for the dimer, Eq. 65. As any state, the ground state
is $2^{N/2}=2$ fold degenerate and identified in this case by $p^{\alpha}=1$
or $p^{\alpha}=-1$ for all $\alpha$. Now consider the effect of a small
perturbation, $H_{N=2}\rightarrow H_{N=2}+vp^{x}$. This does not correspond to
any physical perturbation in terms of spin operators but lifts the ground
state degeneracy. From this point of view, the ground state expectation value
$\expectationvalue{p^{\alpha}}=0$ is fragile, any finite perturbation
violating the conservation of $\tau_{i}$ as defined in Eq. 6 with $i=0,1$
generically causes $\expectationvalue{p^{\alpha}}=\pm 1$. This effect is of
course alleviated at finite temperature, where the relative population
difference of the two lowest states split by $\sim v$ is controlled by the
ratio $v/T$. Kubo’s formula allows to formalize the above considerations for
the linear response of $\expectationvalue{p^{\alpha}}$ with respect to
$vp^{x}$,

$\expectationvalue{p^{\alpha}}=-vG^{R}_{p^{\alpha}p^{x}}(i\omega_{k}=0)\text{.}$
(66)

In Matsubara frequency space, the retarded Green’s function above may be
obtained in the Lehmann representation noting that the parities are diagonal
in the eigenbasis of the unperturbed Hamiltonian
$\matrixelement{n}{p^{\alpha}}{m}=p^{\alpha}_{n}\delta_{nm}$,

$\
G_{p^{\alpha}p^{x}}(i\omega_{k}=0)=\frac{\beta}{\mathcal{Z}}\sum_{n}e^{-\beta
E_{n}}p^{\alpha}_{n}p^{x}_{n}\text{.}$ (67)

At low temperatures this yields $\beta=\frac{1}{T}$, similar to the Curie-like
$1/T$ behaviour of the spin susceptibility of a free spin 1/2 which also
features a degenerate ground state in the field-free case. In complete analogy
to the spin susceptibility in Eq. (63), we can now find the the tree expansion
of the parity susceptibility $G_{p^{\alpha}p^{x}}(i\omega_{k}=0)$ in terms of
the non-local vertices of type $\Gamma_{a}$ (for $\alpha=x$) or $\Gamma_{c}$
($\alpha=y,z$). The expressions are similar to Eq. (63) but crucially probe
different frequency combinations of the vertices ($t=0$ instead of $s=0$). In
other words, non-local vertex components of order $\sim 1/T$ are inherently
expected in the SO(3) Majorana representation. In an exact calculation, these
components are responsible for the $1/T$ parity susceptibility of Eq. (67),
but do not affect the spin susceptiblity. However, the PMFRG is not an exact
method and the unphysical behavior of $\chi^{\Lambda}_{ij}$ at $T=0$ must be a
consequence of truncating the PMFRG flow equations which apparently causes
this divergence to spill over to the spin susceptibility. It is an interesting
question if an improved two-loop truncation scheme (correct to order
$\mathcal{O}(J^{3}))$ [46] or a recently developed but numerically demanding
multi-loop generalizations of the (PF)FRG [62, 63], can be a possible cure to
this problem.

Fortunately, as the unphysical divergence in the PMFRG flow only occurs at
$\Lambda=0$ and for $T=0$, there are other options to extract physically
meaningful results without going beyond the flow equations presented above.
First, it is still possible to detect magnetic phases, heralded by divergences
at finite $\Lambda$ as we have tested for the $J_{1}-J_{2}$ square lattice
Heisenberg model (data not shown).

We devote the rest of the discussion to a second option, which is the
restriction to finite temperatures. As explained above, this can be expected
to suppresses the unphysical divergence and we indeed find all vertices and
flowing susceptibilities converge towards $\Lambda\rightarrow 0$, see lower
inset of Fig. 3 for $T=0.1$.

Figure 3: (Free) energy, heat capacity per spin and static susceptibilities of
the Heisenberg dimer with $J=1$ obtained via PMFRG (red symbols) at
$\Lambda=0$ as a function of temperature. Displayed in solid (dashed) grey
lines are the results obtained by (pseudo-fermion) exact diagonalization, as
well as the finite temperature spin susceptibilities of the PFFRG in black
symbols. Each data point corresponds to a fully converged flow with respect to
$\Lambda$ as demonstrated in the exemplary plot at $T=0.1$ (cf. Fig. 2).

###  VIII.2 Dimer and hexamer at finite temperature

Results for the physical finite-$T$ susceptibility of the dimer at $\Lambda=0$
are shown in Fig. 3. For $T\gtrsim 0.2$, we find a very close agreement
between the susceptibility obtained via PMFRG and the exact result (solid
lines) from Eq. 64. The difference between the exact result and the PMFRG
increases with decreasing temperature, in agreement with the discussion in the
previous subsection. We also show analogous results of the PFFRG, where the
presence of unphysical states seriously compromises the accuracy of the
results at any finite temperature scale. To support this interpretation, we
have also included the results of an exact diagonalization scheme of the
pseudo-fermionic Hamiltonian without projecting out unphysical states, further
referred to as PFED. The close agreement between PFFRG and PFED demonstrates
the problematic impact of unphysical states at finite temperatures which so
far has no known resolution. One approach, the Popov-Fedotov projection
scheme, suppresses unphysical states in exact calculations of observables upon
the introduction of an imaginary chemical potential. However, producing a
quarter-period shift of Matsubara frequencies [69, 57], this option has so far
not been integrated in the PFFRG in a satisfactory manner.

Besides the magnetic susceptibility, our solution of the free energy flow
equation enables us to compute a variety of related thermodynamic observables,
such as the energy per spin and the heat capacity, also displayed in Fig. 3.
We observe good agreement at large enough temperatures. At intermediate scales
$T\simeq 0.5$, the quality of the thermodynamic quantities from the PMFRG
decreases as can be seen most clearly from the overestimation of the energy
per spin or the underestimation of the peak in the heat capacity. These
inaccuracies likely stem from the underestimation of the Majorana self-energy
at small frequencies, a known problem in pseudo-fermion FRG approaches to spin
systems of small dimensionality [32].

Analogous results are obtained for larger spin clusters such as the Heisenberg
hexamer, a hexagon of six equivalent spins with nearest and next-nearest
neighbor interactions, $J_{1}=1$ and $J_{2}=0.5$ respectively. As shown in
Fig. 4, the PMFRG results are in good agreement with ED at not too small
temperatures. The susceptibilities are generally more accurate than the
thermodynamic properties. The susceptibility obtained via PFFRG shows large
deviations from ED results at all temperatures. We emphasize again that small
spin clusters are particularly challenging within the FRG framework since its
built-in mean-field limits are generally not expected to describe such systems
accurately. On the other hand, mean-field approaches perform better in higher-
dimensional systems. The FRG is, hence, expected to reach its full potential
for larger or even infinite systems to which we move on in the following
section.

Figure 4: PMFRG results for the Heisenberg hexamer in analogy to Fig. 3. The
corresponding PFFRG and PFED results of the spin susceptibility are included
in the inset.

##  IX Application: Frustrated spin systems in 2D

We now turn to the application of the PMFRG to two-dimensional, frustrated and
translational invariant Heisenberg spin models described by Hamiltonian (1).
We first study the $J_{1}-J_{2}$ Heisenberg model on the square lattice with
the parameter choice $J_{2}=0.5$ (where the system is expected to be non-
magnetic) and then turn to the triangular lattice model with only nearest
neighbor interaction, $J_{2}=0$. We work at finite temperature $T>0$
throughout and directly in the thermodynamic limit. Thus, as a technical
modification from the previous section, we are required to limit the range of
vertices to $|\mathbf{r}_{i}-\mathbf{r}_{j}|\leq L$, measured in units of the
nearest-neighbor distance [22]. Beyond this distance, vertices (and thus
connected Green functions) are set to zero. We take $L\simeq 10$ large enough
such that our results are converged in $L$. We study the same observables as
in the previous section but report the uniform static susceptibility $\chi/N$
instead of $\chi_{ij}$. In contrast to the previous section, we plot these
observables over $\beta=1/T$.

Figure 5: PMFRG results (dots) for the $J_{1}-J_{2}$ square lattice Heisenberg
model at $J_{2}/J_{1}=0.5$. The panels depict the single site contribution to
the interaction correction to the free energy, internal energy, heat capacity
and uniform susceptibility (top to bottom). The HTSE data (dashed line, up to
9th order) is taken from Ref. [82], its 4,5 Padé approximant is shown as a
solid line. The iPEPS result for the ground state energy $E_{0}/N=-0.495$ from
Ref. [83] is indicated as a dotted line. Figure 6: PMFRG results (dots) for
the nearest-neighbor triangular lattice Heisenberg model. The observables
presented are analogous to Fig. 5. The HTSE data (dashed line, up to 12th
order) is taken from Ref. [84], its 6,6 Padé approximant is shown as a solid
line.

Our PMFRG results for the square lattice are shown in Fig. 5 (dots). We
compare to the high-temperature series expansion (HTSE, dashed line) [82] and
its 4,5 Padé approximant (solid line) with an extended range of stability
$\beta\lesssim 2$ [85], to which our data is in reasonable agreement. We are
not aware of $T>0$ tensor network results for the chosen model, but depict the
iPEPS ground state energy $E_{0}/N=-0.495$ from Ref. [83] (dotted line).
Finally, we remark that when applied to the unfrustrated nearest-neighbor
Heisenberg model ($J_{2}=0$, data not shown), the PMFRG results agree only to
the first order HTSE but deviate strongly from higher order and Monte Carlo
data already for $T=1$. The likely reason is that for the current level of
truncation of flow equations, the FRG is known to violate the Mermin-Wagner
theorem [46], and does, hence, not accurately capture the onset of magnetic
order at $T=0$ in an unfrustrated Heisenberg system.

In Fig. 6, we show the PMFRG results for the triangular lattice nearest
neighbor Heisenberg model (dots). Agreement to the HTSE data [84] (dashed
line, up to 12th order) and its 6,6 Padé approximant is similar as in the
$J_{1}-J_{2}$ square lattice Heisenberg model of Fig. 5. In the temperature
range for which the Padé-HTSE is shown, its accuracy was confirmed by recent
experiments [86] and tensor network results [87].

##  X Conclusion and outlook

In this work, we proposed a FRG approach to spin-1/2 quantum magnets with spin
operators rewritten in the $SO(3)$ Majorana representation. Compared to the
established PFFRG based on representing spins by complex fermions, our PMFRG
method comes with a number of important conceptual differences, both on a
technical level as well as regarding the scope for applications. First, as the
Majorana nature of the spin representation is essential, we derived general
FRG flow equations for generic interacting Majorana Hamiltonians. These could
potentially be useful for other applications [88]. Second, the $SO(3)$
Majorana representation avoids the unphysical states inherent in the complex
fermion representation and instead features a redundant description of spin
states reflected in a fixed artificial degeneracy. As a consequence, the
truncation of flow equations is the only physical approximation made in the
PMFRG. This explains why the PMFRG yields reasonably accurate results for
finite temperatures, being out of reach for the PFFRG. In particular, we
showed how the PMFRG can be used to compute thermodynamic quantities which are
of great experimental relevance. On the downside, the PMFRG’s precision at low
temperatures suffers from a divergence of the $T=0$ flow, which we showed to
be closely related to the (ground-)state degeneracy inherent in the $SO(3)$
Majorana representation, but ultimately caused due to inaccuracies introduced
through the truncation of the hierarchy of flow equations. We thus conclude
that, at the current stage, the PMFRG should be regarded not as a competitor
to the PFFRG, but rather a complement in the practitioners toolbox tailored
for finite and not too small temperatures.

Further work should investigate the potential of the recently proposed
multiloop extension of the (PF)FRG [71, 63, 62] to mitigate the unphysical
divergence mentioned above. Moreover, while the current paper has focused on
Heisenberg systems with global spin rotation symmetry, generalization towards
different classes of systems with reduced symmetries, i.e. Kitaev models and
their variants, should be straightforward. Finally, we emphasize that the
$SO(3)$ Majorana representation is only one out of several Majorana based spin
representations [75]. Based on our results, we believe that these are
promising but relatively underexplored venture points for the application of
many-body methods in the study of spin systems.

## Acknowledgements

We acknowledge useful discussions with Maxime Dupont, Christoph Karrasch,
Alexander Penner, Achim Rosch, and Andreas Weichselbaum. Computations were
performed at the Curta cluster at Freie Universität Berlin and at the
Lawrencium cluster at Lawrence Berkeley National Lab. BS acknowledges
financial support by the German National Academy of Sciences Leopoldina
through grant LPDS 2018-12. NN further acknowledges support from the German
Research Foundation within the TR 183 (project A04).

## References

  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism_ , Springer Series in Solid-State Sciences (Springer, Heidelberg, 2011).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_ , 2nd ed. (Cambridge University, Cambridge, England, 2011).
  * Vojta [2003] M. Vojta, Reports on Progress in Physics 66, 2069 (2003).
  * Anderson [1973] P. Anderson, Mater. Res. Bull. 8, 153 (1973).
  * Balents [2010] L. Balents, Nature (London) 464, 199 (2010).
  * Savary and Balents [2017] L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2017).
  * Wen [1990] X. G. Wen, International Journal of Modern Physics B 04, 239 (1990).
  * Bennett and DiVincenzo [2000] C. H. Bennett and D. P. DiVincenzo, Nature 404, 247 (2000).
  * Sandvik and Kurkijärvi [1991] A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991).
  * Sandvik [1999] A. W. Sandvik, Phys. Rev. B 59, R14157 (1999).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2005] U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Schollwöck [2011] U. Schollwöck, Annals of Physics 326, 96 (2011).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Advances in Physics 57, 143 (2008).
  * Orús [2014] R. Orús, Annals of Physics 349, 117 (2014).
  * Wetterich [1993] C. Wetterich, Physics Letters B (1993), 10.1016/0370-2693(93)90726-X, arXiv:1710.05815 .
  * Kopietz  _et al._ [2010] P. Kopietz, L. Bartosch, and F. Schütz,  _Introduction to the Functional Renormalization Group_ , Lecture Notes in Physics (Springer, Heidelberg, 2010).
  * Polchinski [1984] J. Polchinski, Nuclear Physics B 231, 269 (1984).
  * Metzner  _et al._ [2012] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Rev. Mod. Phys. 84, 299 (2012).
  * Platt  _et al._ [2013] C. Platt, W. Hanke, and R. Thomale, Advances in Physics 62, 453 (2013).
  * Halboth and Metzner [2000] C. J. Halboth and W. Metzner, Phys. Rev. B 61, 7364 (2000).
  * Reuther and Wölfle [2010] J. Reuther and P. Wölfle, Phys. Rev. B 81, 144410 (2010).
  * Abrikosov [1965] A. A. Abrikosov, Physics Physique Fizika 2, 5 (1965).
  * Reuther and Thomale [2011] J. Reuther and R. Thomale, Phys. Rev. B 83, 024402 (2011).
  * Reuther  _et al._ [2011a] J. Reuther, P. Wölfle, R. Darradi, W. Brenig, M. Arlego, and J. Richter, Phys. Rev. B 83, 064416 (2011a).
  * Reuther  _et al._ [2011b] J. Reuther, D. A. Abanin, and R. Thomale, Phys. Rev. B 84, 014417 (2011b).
  * Reuther  _et al._ [2011c] J. Reuther, R. Thomale, and S. Trebst, Phys. Rev. B 84, 100406(R) (2011c).
  * Singh  _et al._ [2012] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Phys. Rev. Lett. 108, 127203 (2012).
  * Reuther  _et al._ [2012] J. Reuther, R. Thomale, and S. Rachel, Phys. Rev. B 86, 155127 (2012).
  * Göttel  _et al._ [2012] S. Göttel, S. Andergassen, C. Honerkamp, D. Schuricht, and S. Wessel, Phys. Rev. B 85, 214406 (2012).
  * Suttner  _et al._ [2014] R. Suttner, C. Platt, J. Reuther, and R. Thomale, Phys. Rev. B 89, 020408(R) (2014).
  * Reuther and Thomale [2014] J. Reuther and R. Thomale, Phys. Rev. B 89, 024412 (2014).
  * Reuther  _et al._ [2014] J. Reuther, R. Thomale, and S. Rachel, Phys. Rev. B 90, 100405(R) (2014).
  * Iqbal  _et al._ [2015] Y. Iqbal, H. O. Jeschke, J. Reuther, R. Valentí, I. I. Mazin, M. Greiter, and R. Thomale, Phys. Rev. B 92, 220404(R) (2015).
  * Rousochatzakis  _et al._ [2015] I. Rousochatzakis, J. Reuther, R. Thomale, S. Rachel, and N. B. Perkins, Phys. Rev. X 5, 041035 (2015).
  * Iqbal  _et al._ [2016a] Y. Iqbal, R. Thomale, F. Parisen Toldin, S. Rachel, and J. Reuther, Phys. Rev. B 94, 140408(R) (2016a).
  * Balz  _et al._ [2016] C. Balz, B. Lake, J. Reuther, H. Luetkens, R. Schönemann, T. Herrmannsdörfer, Y. Singh, A. T. M. Nazmul Islam, E. M. Wheeler, J. A. Rodriguez-Rivera, T. Guidi, G. G. Simeoni, C. Baines, and H. Ryll, Nature Physics 12, 942 (2016).
  * Iqbal  _et al._ [2016b] Y. Iqbal, P. Ghosh, R. Narayanan, B. Kumar, J. Reuther, and R. Thomale, Phys. Rev. B 94, 224403 (2016b).
  * Buessen and Trebst [2016] F. L. Buessen and S. Trebst, Phys. Rev. B 94, 235138 (2016).
  * Hering and Reuther [2017] M. Hering and J. Reuther, Phys. Rev. B 95, 054418 (2017).
  * Baez and Reuther [2017] M. L. Baez and J. Reuther, Phys. Rev. B 96, 045144 (2017).
  * Iqbal  _et al._ [2017] Y. Iqbal, T. Müller, K. Riedl, J. Reuther, S. Rachel, R. Valentí, M. J. P. Gingras, R. Thomale, and H. O. Jeschke, Phys. Rev. Materials 1, 071201(R) (2017).
  * Buessen  _et al._ [2018a] F. L. Buessen, M. Hering, J. Reuther, and S. Trebst, Phys. Rev. Lett. 120, 057201 (2018a).
  * Buessen  _et al._ [2018b] F. L. Buessen, D. Roscher, S. Diehl, and S. Trebst, Phys. Rev. B 97, 064415 (2018b).
  * Roscher  _et al._ [2018] D. Roscher, F. L. Buessen, M. M. Scherer, S. Trebst, and S. Diehl, Phys. Rev. B 97, 064416 (2018).
  * Rück and Reuther [2018] M. Rück and J. Reuther, Physical Review B 97, 144404 (2018).
  * Iqbal  _et al._ [2018] Y. Iqbal, T. Müller, H. O. Jeschke, R. Thomale, and J. Reuther, Phys. Rev. B 98, 064427 (2018).
  * Keleş and Zhao [2018a] A. Keleş and E. Zhao, Phys. Rev. Lett. 120, 187202 (2018a).
  * Keleş and Zhao [2018b] A. Keleş and E. Zhao, Phys. Rev. B 97, 245105 (2018b).
  * Iqbal  _et al._ [2019] Y. Iqbal, T. Müller, P. Ghosh, M. J. P. Gingras, H. O. Jeschke, S. Rachel, J. Reuther, and R. Thomale, Phys. Rev. X 9, 011005 (2019).
  * Hering  _et al._ [2019] M. Hering, J. Sonnenschein, Y. Iqbal, and J. Reuther, Phys. Rev. B 99, 100405(R) (2019).
  * Ghosh  _et al._ [2019a] P. Ghosh, T. Müller, F. P. Toldin, J. Richter, R. Narayanan, R. Thomale, J. Reuther, and Y. Iqbal, Phys. Rev. B 100, 014420 (2019a).
  * Ghosh  _et al._ [2019b] P. Ghosh, Y. Iqbal, T. Müller, R. T. Ponnaganti, R. Thomale, R. Narayanan, J. Reuther, M. J. P. Gingras, and H. O. Jeschke, npj Quantum Materials 4, 63 (2019b).
  * Buessen  _et al._ [2019] F. L. Buessen, V. Noculak, S. Trebst, and J. Reuther, Phys. Rev. B 100, 125164 (2019).
  * Niggemann  _et al._ [2019] N. Niggemann, M. Hering, and J. Reuther, Journal of Physics: Condensed Matter 32, 024001 (2019).
  * Revelli  _et al._ [2019] A. Revelli, C. C. Loo, D. Kiese, P. Becker, T. Fröhlich, T. Lorenz, M. Moretti Sala, G. Monaco, F. L. Buessen, J. Attig, M. Hermanns, S. V. Streltsov, D. I. Khomskii, J. van den Brink, M. Braden, P. H. M. van Loosdrecht, S. Trebst, A. Paramekanti, and M. Grüninger, Phys. Rev. B 100, 085139 (2019).
  * Roscher  _et al._ [2019] D. Roscher, N. Gneist, M. M. Scherer, S. Trebst, and S. Diehl, Phys. Rev. B 100, 125130 (2019).
  * Chillal  _et al._ [2020] S. Chillal, Y. Iqbal, H. O. Jeschke, J. A. Rodriguez-Rivera, R. Bewley, P. Manuel, D. Khalyavin, P. Steffens, R. Thomale, A. T. M. N. Islam, J. Reuther, and B. Lake, Nature Communications 11, 2348 (2020).
  * Revelli  _et al._ [2020] A. Revelli, M. Moretti Sala, G. Monaco, C. Hickey, P. Becker, F. Freund, A. Jesche, P. Gegenwart, T. Eschmann, F. L. Buessen, S. Trebst, P. H. M. van Loosdrecht, J. van den Brink, and M. Grüninger, Phys. Rev. Research 2, 043094 (2020).
  * Kiese  _et al._ [2020a] D. Kiese, F. L. Buessen, C. Hickey, S. Trebst, and M. M. Scherer, Phys. Rev. Research 2, 013370 (2020a).
  * Iida  _et al._ [2020] K. Iida, H. K. Yoshida, A. Nakao, H. O. Jeschke, Y. Iqbal, K. Nakajima, S. Ohira-Kawamura, K. Munakata, Y. Inamura, N. Murai, M. Ishikado, R. Kumai, T. Okada, M. Oda, K. Kakurai, and M. Matsuda, Phys. Rev. B 101, 220408(R) (2020).
  * Kiese  _et al._ [2020b] D. Kiese, T. Müller, Y. Iqbal, R. Thomale, and S. Trebst, arxiv (2020b), arXiv:2011.01269 .
  * Thoenniss  _et al._ [2020] J. Thoenniss, M. K. Ritter, F. B. Kugler, J. Von Delft, and M. Punk, arxiv (2020), arXiv:2011.01268 .
  * Martin [1959] J. L. Martin, Proc. R. Soc. London A 251, 536 (1959).
  * Tsvelik [1992] A. M. Tsvelik, Phys. Rev. Lett. 69, 2142 (1992).
  * Kitaev [2006] A. Kitaev, Annals of Physics 321, 2 (2006).
  * Krieg and Kopietz [2019] J. Krieg and P. Kopietz, Phys. Rev. B 99, 060403(R) (2019).
  * Goll  _et al._ [2019] R. Goll, D. Tarasevych, J. Krieg, and P. Kopietz, Phys. Rev. B 100, 174424 (2019).
  * Popov and Fedotov [1988] V. N. Popov and S. A. Fedotov, Sov. Phys. JETP 67, 535 (1988).
  * Kugler and von Delft [2018a] F. B. Kugler and J. von Delft, Phys. Rev. Lett. 120, 057403 (2018a).
  * Kugler and von Delft [2018b] F. B. Kugler and J. von Delft, Phys. Rev. B 97, 035162 (2018b).
  * Kugler and von Delft [2018c] F. B. Kugler and J. von Delft, New Journal of Physics 20, 123029 (2018c).
  * Shnirman and Makhlin [2003] A. Shnirman and Y. Makhlin, Phys. Rev. Lett. 91, 207204 (2003).
  * Biswas  _et al._ [2011] R. R. Biswas, L. Fu, C. R. Laumann, and S. Sachdev, Phys. Rev. B 83, 245131 (2011).
  * Fu  _et al._ [2018] J. Fu, J. Knolle, and N. B. Perkins, Phys. Rev. B 97, 115142 (2018).
  * Schad  _et al._ [2015] P. Schad, Y. Makhlin, B. Narozhny, G. Schön, and A. Shnirman, Annals of Physics 361, 401 (2015).
  * Nilsson and Bazzanella [2013] J. Nilsson and M. Bazzanella, Phys. Rev. B 88, 045112 (2013).
  * Behrends  _et al._ [2019] J. Behrends, J. H. Bardarson, and B. Béri, Phys. Rev. B 99, 195123 (2019).
  * Katanin [2004] A. A. Katanin, Phys. Rev. B 70, 115109 (2004).
  * Morris [1994] T. R. Morris, Int. J. Mod. Phys. A 09, 2411 (1994).
  * Karrasch  _et al._ [2008] C. Karrasch, R. Hedden, R. Peters, T. Pruschke, K. Schönhammer, and V. Meden, Journal of Physics: Condensed Matter 20, 345205 (2008).
  * Rosner  _et al._ [2003] H. Rosner, R. R. P. Singh, W. H. Zheng, J. Oitmaa, and W. E. Pickett, Phys. Rev. B 67, 014416 (2003).
  * Poilblanc and Mambrini [2017] D. Poilblanc and M. Mambrini, Phys. Rev. B 96, 014414 (2017).
  * Elstner  _et al._ [1993] N. Elstner, R. R. P. Singh, and A. P. Young, Phys. Rev. Lett. 71, 1629 (1993).
  * Poilblanc  _et al._ [2020] D. Poilblanc, M. Mambrini, and F. Alet, arxiv (2020), arXiv:2010.07828 .
  * Cui  _et al._ [2018] Y. Cui, J. Dai, P. Zhou, P. S. Wang, T. R. Li, W. H. Song, J. C. Wang, L. Ma, Z. Zhang, S. Y. Li, G. M. Luke, B. Normand, T. Xiang, and W. Yu, Phys. Rev. Materials 2, 044403 (2018).
  * Chen  _et al._ [2019] L. Chen, D.-W. Qu, H. Li, B.-B. Chen, S.-S. Gong, J. von Delft, A. Weichselbaum, and W. Li, Phys. Rev. B 99, 140404(R) (2019).
  * Rahmani and Franz [2019] A. Rahmani and M. Franz, Reports on Progress in Physics 82, 084501 (2019).

##  Appendix A Details on the numerical implementation

The flow equations presented above can be solved using standard, error
controlled Runge-Kutta schemes, such as the fifth-order Dormand-Prince method.
In our case, we found little dependence of our results on the choice of the
integration method used upon decreasing the relative and absolute accuracy to
$\sim 10^{-2}$ or lower. In equivalence to implementations of the PFFRG, the
maximum distance treated in four-point vertices $\Gamma_{ij}$ is limited to
$|\mathbf{r}_{i}-\mathbf{r}_{j}|\leq L\simeq 10$ for translation invariant
systems.

At finite temperatures, we treat the frequency dependence by generating a set
of $N_{\omega}=32$ positive Matsubara indices such that our results are
converged in $N_{\omega}$. The indices were chosen according to the following
scheme such that the smallest frequencies are included exactly, while larger
indices are more sparse and require for linear interpolation in between them:

$n_{i}=\text{round}\left[z\sinh\left(\frac{i}{z}\right)\right],\quad
i=0,1,\dots,N_{\omega}\text{.}$ (68)

The parameter $z$ is then fully determined upon specification of the
temperature, the number of frequencies, and the maximum frequency. Since the
according (fermionic) frequencies are given by $\omega_{n}=\pi T(2n+1)$, one
needs to be careful when implementing fermionic symmetries such as
$\gamma(-n_{\omega})=-\gamma(n_{\omega}-1)$. Furthermore, the Matsubara
integers corresponding to sums and differences of fermionic frequencies are

$\displaystyle\omega_{1}+\omega_{2}$ $\displaystyle\leftrightarrow
n_{\omega_{1}}+n_{\omega_{2}}+1$ $\displaystyle\omega_{1}-\omega_{2}$
$\displaystyle\leftrightarrow n_{\omega_{1}}-n_{\omega_{2}}\text{.}$ (69)

As a result, only those sets of Matsubara integers that sum up to odd integers
$n_{s}+n_{t}+n_{u}=2n_{\omega_{1}}+1$ are physical within energy conservation
and will be evaluated in vertices. For the less robust implementation at
$T=0$, we choose a logarithmic frequency mesh consisting of $N_{\omega}=96$
positive frequencies to avoid numerical errors from the finite frequency grid.
The frequency integral in the Katanin contribution is then carried out
numerically using a trapezoidal method.