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

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 renders 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, particularly, 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.

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 princi-ple, 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, and has constantly been extended and generalized. Today, the PFFRG is, hence, remarkably flexible with a scope of applicability comprising two dimensional [22, 24-35, 37, 38, 40, 41, 44-46, 48, 49, 51, 54, 55, 57, 59-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 arXiv:2012.14836v4 [cond-mat.str-el] 20 Sep 2021 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 will be considered, where (i, j) refers to all possible pairings of sites and S α i is the α 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 =↑, ↓) [23]: Here and in the following, we set = k B = 1. However, this representation is a valid rewriting of the spin operators only in the local subspace with a f † ia 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 PF-FRG may be faithfully implemented with the simpler condition a f † ia f ia = 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 Λ, introduced in the bare Green function G 0,Λ (iω) = G 0 (iω)Θ(|ω|−Λ), is varied. The basic idea is that at the starting point Λ = ∞, the bare propagator vanishes and all vertex functions are trivially known. For a numerical solution of the flow equations down to Λ = 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 Λ 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) → SU (N )] and the spin length [1/2 → S], respectively. In two very recent works, certain contributions of the six-point vertex have been taken into account using a multiloop extension [62,63] 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 α i at site i, three different flavors α ∈ {x, y, z} of Majorana fermions η α † i = η α i are introduced. They fulfill the anticommutation relations {η α i , η β j } = δ ij δ αβ which imply (η α i ) 2 = 1/2. The formal Hilbert space dimension per Majorana is √ 2 as appropriate for half a (complex) fermion. The spin operators S α i = − i 2 βγ ε αβγ η β i η γ i , more explicitly written as can be easily checked to fulfill the spin-1/2 algebra As an example, a Heisenberg coupling term from Hamiltonian (1) is represented as As usual for auxiliary particle representations, the SO(3) Majorana representation comes with a gauge freedom. The local Z 2 gauge transformation η α i → ε i η α i with ε i = ±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 which anticommutes with any τ j from a different site j = i and fulfills Consequently, τ i commutes with all spin operators and thus with any spin Hamiltonian. To construct a set of mutually commuting operators one needs to pair τ i with another conserved Majorana operator. One choice [73] is to define an additional Majorana η 0 i per site, so that the parity p i = 2iτ i η 0 i with eigenvalues ±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-toone correspondence to the original local spin Hilbert space. To invoke η 0 i 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τ i τ j = ±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 Z 2 gauge symmetry: As the transformation τ i → −τ 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 Z pm = 2 N/2 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 where the first term f pm ≡ − T N log (Z pm ) 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) O s are easily computed in the Majorana representation as well. This follows from the observation that the Majorana version of such an operator, O pm , is diagonal in the parity sector and the same is true for any physical density matrix ρ pm , like for example the Boltzmann factor ρ pm ∼ e −βHpm . Then the degeneracy factor 2 N/2 simply cancels [76] and we have 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 τ i , the spin operators can be re-expressed as We may now consider the general SO(3) transformation η α i → β R αβ η β i with R αβ ∈ SO(3) being a three dimensional rotation matrix. As τ i is invariant under this transformation [75], spin operators must transform as It follows that physical SO(3) rotations of a spin i are equivalent to rotations of the Majorana vector (η x i , η y i , η z i ).

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 where {µ i } is an arbitrary set of single-particle indices. As above, we use the convention η µi , η µj = δ µiµ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 η µ (τ ). The action reads where ∂ τ denotes a derivative with respect to imaginary time and β = 1/T . We define the Fourier transform η µ (τ ) = T n e iωnτ η µ (iω n ) where the fermionic Matsubara frequencies are given by iω n = πT (2n + 1), with n ∈ Z. In slight abuse of notation, in the following, we will denote ω n1 by ω 1 and equivalently for other frequencies. The non-interacting part of the action may then be written as (14) with the bare Majorana Green's function 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η 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 Ψ containing an arbitrary number of bosonic or Grassmann fields labeled by the composite index l = (ω l , µ l ), where l = β −1 ω l µ l . A comparison of Eq. (16) and Eq. (14) yields the direct correspondence Ψ l=(µ l ,ω l ) = η µ l (ω l ). We emphasize the difference to the superfield vectors of complex fermions or bosons, which require an additional but independent superfield label, i.e. Ψ = (ψ, ψ).
The starting point of the FRG scheme is the introduction of a cutoff scale Λ in the bare Green's function G 0 → G Λ 0 such that G Λ=∞ 0 = 0 and G Λ=0 0 = G 0 . Although the flow equations describing the evolution of irreducible vertices with Λ [17] below are general, in the rest of this work, we will consider a multiplicative Matsubara frequency cutoff Θ Λ (ω 1 ) to the bare Green's function At zero temperature, this cutoff is often chosen to be a Heaviside function Θ Λ (|ω|) = θ(|ω| − Λ), 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 Σ and the four-point vertex Γ acquire implicit dependence on Λ. These quantities are defined via the Dyson equation in a superspace spanned by (ω i , µ i ) and the tree-expansion for the connected Green's functions respectively. This Λ-dependence is given by coupled differential equations, referred to as flow equations. Physical results can be extracted from the solution at Λ = 0. Since the action for Majorana systems was rephrased in superfield notation, we can employ the associated general flow equations [17] for Σ Λ and Γ Λ . 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, With the above definition, the Dyson equation for fixed fre- , can be written as The Green's function and self-energy defined in Eq. (20a) and 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 The Majorana flow equations for the interaction correction to the free energy, self energy and the four-point vertex read [17] d dΛ 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 In order to solve the flow equations, initial conditions for selfenergy and the four-point vertex are required. As the bare propagator vanishes in this limit, we immediately see that

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), (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 = G 4 µ1,µ2,µ3,µ4 (s, t, u)βδ ω1+ω2+ω3+ω4,0 .
where the labels (1, 2, 3, 4) contain all arguments that are not explicitly specified, i.e 1 = (µ 1 , ω 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 Σ and Γ due to their relations via Eqs. (18) and (19).

A. Hermiticity
The Hamiltonian is a hermitian operator, satisfying H = H † . Due to O * = O † and η(τ ) † = e −Hτ ηe Hτ = η(−τ ) in the Heisenberg picture, one can find the complex conjugate of the two-point Green's functions as 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

B. Time reversal symmetry
Time reversal T is an anti-unitary operation ( ψ|ψ * = T ψ|T ψ ) which in the present case can be defined by performing a complex conjugation while leaving Majorana operators invariant [78]: This flips the sign of the spin operators (3) Similarly, the four-point correlator has the property Since our considerations from here on require the explicit specification of site indices, we will now separate the previously used superlabel µ into a site-index and a Majorana flavor µ → (i, α). In the SO(3) Majorana representation spins are invariant under the gauge transformation η α i → ε i η α i for all α = x, y, z with ε i = ±1 for an arbitrary lattice site i.
Since expectation values must be invariant under gauge transformations as well, we may write where ε i1 ε i2 = −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 Likewise, the four-point correlator can only depend on up to two distinct sites only, so we choose 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 Γ 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].

D. 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 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 Γ 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, Table I. Symmetry transformations corresponding to specific spin rotations along the x, y and z axes.
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.

E. 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 × 3 rotation matrix R αµ (φ) act on the Majorana fermions as 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 I. The combination R x (π/2)•R z (π/2) ≡ P realizes an anti-cyclic permutation of the flavors. We apply these symmetries to correlators, using the convention γ = α = β = γ 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.
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 π/2 transform different flavor combinations into each other, for instance 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 η x Majoranas then transform as η x i → cos θη x i − sin θη y i so that (41) Expanding the product and using the above symmetries, we obtain a relation independent of θ 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 P and rotations around an arbitrary axis may always be decomposed as a product of x, y and z rotations.
In the special case i = j, there are only two independent vertices since Vertices with negative bosonic frequencies are symmetry related to positive frequencies by time-reversal and a symmetry t ↔ u further allows to reduce the numerical effort. Details are given in Table II. In the above parametrization, the flow equations for the interaction correction to the free energy per Operation Symmetry for Γµ ij (s, t, u) valid µ c   Table II. Transformations of the frequency arguments under time reversal T and specific permutations of indices inΓ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 ↔ u symmetry for Γc can be established by the exchange 1 ↔ 2, which changes the vertex to the form Γxyyx. Using Eq. (42) to express Γxyyx(s, t, u) = −Γc(s, u, t) in terms of the other vertices used in the parametrization, we obtain Γc ij (s, u, t) = (−Γa ij + Γ b ij + Γc ij )(s, t, u).
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 Similarly, we may now express the flow equations for fourpoint 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, To reduce the length of expressions, we have defined the single-channel contributions X Λ a,b,c ij andX Λ a,b,c,d ij in Eqs. (51) and (52) [46]. The flow equations of local vertices are obtained noting thatX Λ a,b,c ii (s, t, u) ≡ X Λ a,b,c ii (s, t, u). We further stress that no flow equation for Γ c ii is required in Eq. (50), as this vertex is equivalent to Γ 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 At zero temperature, frequencies become continuous and T ω → (2π) −1 dω. Using the sharp frequency cutoff G 0Λ (ω) = G 0 (ω)θ(|ω| − Λ), we thus obtain in the usual way using Morris's Lemma [80] At finite temperatures, a sharp cutoff of frequencies is no longer possible due to ambiguities that arise if |ω| − Λ 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 Using Eqs. (17), (18) and (53) the expressions for the Green's function and the single-scale propagator become 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 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 Λ J down to Λ 0, approximating the initial conditions with the Λ → ∞ 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 Λ = 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 PM-FRG flow equations.
From the partition function of a N -spin system with eigenenergies E n , Z = n e −βEn , the free energy per spin is given by The energy per spin is which as a function of T also determines the heat capacity 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 Λ=0 int . The non-interacting free energy for three pseudo-Majoranas per site is f pm,0 = −3T log(2)/2. Using the relation between f pm and f , Eq. (8), we finally obtain The static spin-spin correlator can be computed from Note that χ 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 Λ, Of particular interest for the two-dimensional systems below is the uniform susceptibility χ = i,j χ ij .

VIII. APPLICATION: SMALL SPIN CLUSTERS
A. 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 nontrivial 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 = α S α 0 S α 1 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ν + E n − E m = 0 is often neglected in textbook derivations. We obtain 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 χ Λ ij flows smoothly without any feature, surpasses the exact results χ ij = ±0.5 and diverges at Λ = 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 Γ a,01 = Γ x0,x0,x1,x1 and Γ c,01 = Γ 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, Here, p α ≡ 2iη α 0 η α 1 are the three flavor parities related to the non-local parity introduced in Sec. III via p (0,1) = 2iτ 0 τ 1 = −p x p y p z . While p (i,j) = ±1 is always conserved for generic spin systems, p α = ±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 α = 1 or p α = −1 for all α. Now consider the effect of a small perturbation, H N =2 → 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 p α = 0 is fragile, any finite perturbation violating the conservation of τ i as defined in Eq. (6) with i = 0, 1 generically causes p α = ±1. This effect is of course alleviated at finite temperature, where the relative population difference of the two lowest states split by ∼ v is controlled by the ratio v/T . Kubo's formula allows to formalize the above considerations for the linear response of p α with respect to vp x , 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 n|p α |m = p α n δ nm , At low temperatures this yields β = 1 T , similar to the Curielike 1/T behaviour of the spin susceptibility of a free spin 1/2 which also features a degenerate ground state in the fieldfree 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 α p x (iω k = 0) in terms of the non-local vertices of type Γ a (for α = x) or Γ c (α = 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 ∼ 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 χ Λ 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 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 Λ = 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 Λ 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 Λ → 0, see lower inset of Fig. 3 for T = 0.1.

B. Dimer and hexamer at finite temperature
Results for the physical finite-T susceptibility of the dimer at Λ = 0 are shown in Fig. 3. For T 0.2, we find a very close agreement between the susceptibility obtained via PM-FRG 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 [57,69], 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 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 pseudofermion 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.

IX. APPLICATION: FRUSTRATED SPIN SYSTEMS IN 2D
We now turn to the application of the PMFRG to twodimensional, 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 |r i − r j | ≤ 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 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 χ/N instead of χ ij . In contrast to the previous section, we plot these observables over β = 1/T . 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 β 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 PF-FRG 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 [62,63,71] 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 straightfor-ward. 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.