Reference results for the momentum space functional renormalization group

The functional renormalization group (FRG), an established computational method for quantum many-body phenomena, has been subject to a diversification in topical applications, analytic approximations and numerical implementations. Despite significant efforts to accomplish a coherent standard through benchmarks and the reproduction of previous results, no systematic and comprehensive comparison has been provided until now. While this has not prevented the publication of relevant scientific results we argue that established mutual agreement across realizations will strengthen confidence in the method. To this end, we report explicit implementational details and numerical data reproduced thrice independently up to machine accuracy. To substantiate the reproducibility of our calculations, we scrutinize pillar FRG results reported in the literature, and discuss our calculations of these reference systems. We mean to entice other groups to reproduce and establish this set of benchmark FRG results thus propagating the joint effort of the FRG community to engage in a shared knowledge repository as a reference standard for FRG implementations


Introduction
The functional renormalization group (FRG) is one of the most promising methods for the analysis of lowtemperature instabilities of two dimensional materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] as it offers an unbiased view, not presupposing the existence of certain phase transitions [17].Because the numerics involved in prediction using the FRG are expensive, the full breadth of possible dependencies could only be captured by a few implementations of the twodimensional single-band Hubbard model on the square lattice [18,19], while a generalization to non-SU (2) symmetric or multi-orbital models is still absent.This inherent fragmentation has rendered it difficult to directly compare distinct implementations.Traditionally only qualitative results have been employed [20].Unfortunately various qualitative features produced by the FRG, such as the critical scale and leading instabilities, have (over the course of converging our three distinct implementations) proven to be robust against a certain degree ‡ These authors contributed equally a E-mail: beyer@physik.rwth-aachen.deb E-mail: Profe@itp.uni-frankfurt.dec E-mail: klebl@physik.rwth-aachen.de of programming error e.g. the exact ordering of indices in subleading diagrams.When trying to predict manybody instabilities of novel systems, certainty over reference calculations is crucial to assure validity.Beyond the qualitative level, the aspiration of the method is to become a quantitative quantum many-body tool, further emphasizing numerical correctness.The definition of correctness here refers to a certain formal approximation level, i.e. a specific truncation of the hierarchy of FRG equations, as detailed in Ref. [17].When gauging the physical potential of a certain approximation or comparing various approximations, validity of numerical implementations is mandatory.It is therefore desirable to introduce a shared knowledge of the "proper" results of the tensor contractions intrinsic to most FRG calculation.A reproduction of the exact numerics yields certainty over prefactors, signs as well as contraction indices, not ensured when validating using e.g. the phase transitions predicted in the square lattice Hubbard model.We try to lay the groundwork for achieving consistency across different codes by providing momentum-space data of the frequency independent interaction vertex reproduced by three different, independently developed FRG codes.We choose this focus as a first step towards a shared knowledge repository, motivated by the application to strongly correlated states in two-dimensional materials.
The paper is structured as follows: We first give a short briefing on the theory and some technicalities of the momentum-space FRG and its truncated unity formfactor expansion.Then we introduce the three distinct implementations and elaborate their areas of application.We proceed to provide results in the three tests systems we have chosen, in order to establish agreement with published references.The penultimate section presents the exact procedure of our comparison, enabling the reader to compare their own implementation using the provided parameter and data sets.Intricacies and specifics regarding various challenges faced during development are presented in the Appendix.

Functional Renormalization Group
The derivation of FRG is given in utmost brevity here, we refer the interested reader to Refs.[17,21] for an in-depth discussion.We commence with the partition function of the interacting system: where ψ, ψ are fermionic Grassmann fields and S is the action of the system given by single-particle and interacting components: Here Q 0 = ik 0 − H 0 + µ = G −1 0 denotes the inverse of the bare one-particle Green's function with H 0 as the non-interacting part of the Hamiltonian, k 0 a Matsubara frequency and µ the chemical potential.
We can now take a derivative of the effective action with respect to the scale introduced by means of the regulator and obtain the Wetterich equation [22,23]: with and Expanding in the order of the fields ψ, ψ we obtain an infinite hierarchy of differential equations with Γ Λ,(n) dependent on all even contributions from Γ Λ, (2) to Γ Λ,(n+2) .Next let us clarify the truncation level used for the benchmarks presented in this publication: We note that higher-order derivatives become decreasingly relevant because they correspond to multi-electron interactions [24].Thus, we truncate at Γ Λ, (4) , neglecting all higher order contributions of the Grassmann field expansion.Furthermore we neglect the self-energy contribution Γ Λ, (2)  as -within the scope of this work -Matsubara frequency dependencies of the vertex functions are ignored [25].
Their inclusion is possible and has been covered extensively for the 2D Hubbard model in Refs.[18,19,[26][27][28][29], considering the comparison sought here they are however inopportune.We express these equations as diagrams with Γ Λ,(4) = V 1 as 4-particle nodes connected by the G 0 propagators and their derivative: The remaining equation is then diagrammatically given in Fig. 1.It is insightful to note that up to the truncation of Γ (2) and Γ (n≥6) in the Grassmann field expansion the FRG introduces no approximations.
To obtain a prediction for the two-particle interaction V independent of the artificial scale we integrate the ordinary differential flow equation (Fig. 1) starting at infinite (large compared to bandwidth) Λ and iterating towards Λ = 0.When encountering a phasetransition, the equations will diverge [24,30], making Γ Λ, (6) a meaningful contribution and invalidating the above truncation.We terminate the integration at this point and examine the final V Λc to determine the ordering associated to the phase transition.While this approximation can be improved upon via the inclusion of more diagrams, e.g. by applying the more elaborate multiloop-FRG or Parquet-approximations [18,19,[31][32][33], for ease of comparison and reduction in computational cost it is advantageous to maintain ease of access.

Regulators
For the calculations in this work two different types of regulators are used.The sharp-frequency-cutoff Θ Λ (k) = θ(|k 0 | − Λ) and the Ω-cutoff [34]: ).Since we assume a static vertex and neglect selfenergy effects we can analytically calculate all Matsubara frequency summations needed during the integrations depending on the type of regulator chosen.Through this we obtain the following loop derivatives for the sharp-cutoff: where we used the abbreviated notation = b1 (l) and = b2 (l ) (b 1 and b 2 refer to band indices).In the Ω-cutoff scheme, using the same notation the resulting equations for the loop derivative are given by: LΛ,b1b2 In either case the upper (+) sign corresponds to the particle-hole loop in the C-and D-channel diagrams while the lower (−) sign corresponds to the particleparticle loop in the P -channel diagram.In the case of the Ω regulator we neglect the evaluation of the additional poles of the Matsubara summation which might occur when Λ = , Λ = or Λ = = .This is justified by calculation of the dispersion only on a discrete grid, which delegates these cases to be extremely unlikely.

Band vs Orbital Calculation
The above mentioned diagrams (cf.Fig. 1) can be evaluated in either band or orbital space.To illustrate the discrepancies between the two we here revert to the nondiagrammatic representation of the P -channel.In this channel we can specify l to be: We must calculate the loop derivative in band space yielding Lb1b2 − (l, l ).If we desire to integrate the momentum dependence l using objects represented in orbital space we need to transform this according to (9) in each flow iteration.Alternatively we can transform the initial vertex into band space and revert to orbital space using and its inverse after the FRG flow.The resulting flow equations in band or orbital space respectively are given as For convenience we have introduced Fermi surface single patch Fig. 2 Meshing schemes for the functional renormalization group calculations.(a) N -patch meshing scheme (figure adapted from Ref. [35]), the points are defined at the intersections of the Fermi surface (red) with the patch-lines (dashed gray).For a single patch, the yellow, filled area is projected to the same momentum point.(b) Bravais momentum mesh in primitive zone (PZ) instead of Brillouin zone (BZ) with mesh refinement around one selected point shown in orange.
the properly normalized integrals over Brillouin zone (BZ) and primitive zone (PZ).The computational cost can be reduced by using band space representations of both vertex and loop, but because Γ Λ, (4) is not an observable, this might introduce gauge phases into the calculations.These can lead to difficulties when interpreting the results and are absent when calculating in orbital space.The calculations within this work have therefore been done in orbital representation, where applicable.Within the scope of FRG we treat all non-diagonal quantum numbers on equal footing and use one shared, linearized index o i to describe spin, orbitals, sites, etc.Having dealt with the discrete quantum numbers we now turn our attention to the more complicated, artificially discretized representation of the continuous momentum dependencies of the vertex k.

Grid-FRG and Refinement
Since the following section will deal with different methods for the evaluation of the momentum dependencies we start with an understanding of the fundamental requirements posed and introduce the basic methods of patching the BZ.
The momentum resolution of the functional renormalization group is limited by two distinct issues: The lower bound is given by the ability to resolve sharp features of the fermionic momenta in the loop derivative at low scale.As these occur in the proximity of the Fermisurface (FS) this is equivalent to demanding an accurate representation of said FS within the chosen momentum scheme.The upper bound is determined by computational limitations.Because the memory requirement for saving the vertices -the largest objects -scales pro-portional to N 3 k , raising the fidelity of momentum discretization is constrained by availability of RAM.
The traditional method of fulfilling these two aspects is to allocate a number of points along the FS of the system to be used as momentum points for the vertex calculations [see Fig. 2 (a)].While the vertex would only be defined on the FS, using the assumption that it is almost constant within the patch [35][36][37], the propagators and their derivatives are evaluated for a set of points spanning the entire area of the patch.The values would be averaged to obtain the value of L at the FS patch point, finely sampling the patches area.This method is discussed in more detail in Ref. [35].
The immediate issue with this parameterization is the evaluation of the fourth vertex momentum.Having constrained k 1 , k 2 , k 3 to points on the FS we have no capacity to guarantee that k 4 = k 1 + k 2 − k 3 coincides with one of the patch centers.The assumption is therefore made that this can be projected to the closest existing point, the center of the patch we land within.This method, while historically successful in the prediction of fundamental models [8,10] breaks momentum conservation of the interaction.This is especially problematic when attempting to analyze spin-orbit-coupled (SOC) or multi-site systems where conserving momentum becomes important due to momentum-spin-locking or momentum-orbit-locking.
The easiest method to restore momentum conservation is the introduction of an equispaced ("Bravais") mesh within the first PZ.For this case, shown in Fig. 2 (b), the addition/subtraction of any number of momentum points is assured to be a valid momentum point.The obvious issue though is that the majority of points is now far from the FS, making the resolution of the aforementioned features doubtful.
If we were able to use arbitrarily high momentum resolutions this method would yield the exact momentum integrations.It therefore is the base-case for the simulations performed in this paper and offers the easiest implementation of all two dimensional momentum space FRG alternatives.The main challenges to overcome are the memory constrains imposed by the scaling with momentum resolution and orbital degrees of freedom.To circumvent this issue we introduce an averaging similar to the one performed in the evaluation of N -patch as where l is taken from the normal (coarse) fermionic mesh and l f is a small (refined) offset relative to this point.We choose the sign to conserve total momentum in each channel.The resulting refined grid is shown in Fig. 2 (b).Crucially, we employ the refinement only for an averaging within the calculations of the loop derivative and require no information of the vertex on the full set of refined momentum points.As in N -patch, we assume the interactions to be sufficiently slowly varying such that V (k This approximation is even more reasonable here as the maximum distance from the center is greatly decreased.Some intricacies on the exact implementation of choosing i, j, m, n from orbital or band space are discussed in Appendix E and are omitted here.

Truncated Unities
Building upon the presented basic implementation of grid-FRG we now expand to further approximations.The truncated unity extension to the functional renormalization group (TUFRG) is an approach to reduce the computational complexity introduced using grid-FRG while maintaining momentum conservation.Noticing the divergences introduced in the FRG flow will stem from the loop derivative, we isolate the bosonic momentum q of L± (l, ±(q − l)).For each channel (P , C, D) in Fig. 1 we thus obtain a distinct transfer momentum: The remaining weaker fermionic dependencies of vertex and loop are expanded in a (symmetrized) basis of linearly independent functions f m (k) -called formfactors.The choice of these functions is based on the construction of basis functions for irreducible representations described in Ref. [38] with a few improvements and generalizations detailed in Appendix G.We project/unproject the flow equations given in Fig. 1 into/from the formfactor basis using the following operators for each channel X ∈ P , C, D: It is instructive to note that with a complete basis of formfactors this is an exact unitary transformation, but in our case we truncate the basis by neglecting longrange formfactors.This is reasonable because the interaction strength declines at long ranges [39].The truncation makes the truncated unity expansion an approximation to the grid-FRG implementation, but significantly reduces the memory requirements for the vertices which shrink by a factor of (N f /N k ) 2 .This method has been developed in Ref. [39] based on the earlier work in Ref. [9] and has proven successful for Hubbard and graphene type models [1,3,4,39].When considering the vertex in band space, the gauge phases carried by the orbital-to-band matrices can not be disregarded.Thus, the band space vertex cannot be captured accurately in the truncated unity projections, this approximation thus necessitates orbital space calculations.TUFRG also has slight symmetry breaking at long coupling ranges for multi-site models, the details and workaround of this are described in Appendix F.

Real Space Truncated Unities
In models with broken translational symmetry, momentum is not good quantum number and we thus require a different approach to treat these systems.The proposal here is to use the real space equivalents of formfactors, which we will refer to as bonds for the sake of clarity.In real space the bosonic momentum dependence of each channel is translated to a dependence on two orbitals (or sites).The fermionic momenta are then equivalent to the remaining orbital dependencies.We define the real space bonds g bi (r j ) as Kronecker deltas spanning the lattice where r i is the position of the site i and b i is a connection to another site starting at site i.In the case of large unit cells, where momentum is still conserved, it may also be beneficial to introduce real space truncated unities.The projections in this mixed space representation then follow Eqs.(18a) to (18c): The generation of the sets of real space bonds and momentum space formfactors can now be performed following one of two different rules: For the first, simpler one, we assume that the two sets have no connection, allowing for a two step projection with separate real-and momentum space unities.This decouples the implementations but comes with a severe drawback: For models with more than a single site per unit cell, it is impossible to include the same interaction orders for each site, thus breaking the rotational and inversion symmetries as detailed in Appendix F. The second approach remedies these concerns by first creating all real space bonds on the full lattice with n l1 •n k1 ×n l2 •n k2 elements, where n li is the number of real space unit cells and n ki is the number of momentum points in the direction of basis vector i.Afterwards we perform a Fourier transform of all bonds and obtain formfactor-bonds of the following general form where B bi is the beyond unit cell part of the bond b i .
The real space TU was first developed for one dimensional systems starting from an analysis of the perturbative orders [40][41][42] and was thereafter re-derived for arbitrary dimension in the context of truncated unities as we present it here [43].

Implementations
Having explained the theoretical backdrop for the implementations we now want to discuss the three distinct FRG codes independently developed by the authors, sharing minor details of the code-base (input routines, output routines).Each of the programs will be given a short introduction to establish the functionalities it covers, its merits and the contrast to the other implementations.

Code #1: grid-FRG
As the generation of band structures from ab initio simulations is common practice, it is advantageous to be able to use these as starting point for FRG calculations.Thus, this code implements the static four-point FRG equations in either band or orbital space and on a regular momentum grid and allows the study of correlation effects of arbitrary periodic systems if the following can be provided: -A few isolated low-energy bands of the material in momentum space on a regular, fine momentum mesh.-Addition and subtraction rules for the momentum mesh.-The 4-point vertex in any basis that can be connected to the band basis by unitary transformation (and thus is of the same dimension).As V may be an object of large size, it is possible to supply it on a coarse momentum mesh.-The orbital to band transform [Bloch functions u ob (k)] for all momentum points.In case the 4-point vertex is given in band space, the matrices fulfill û(k) ≡ 1.
Since the code is designed to operate on very general models defined only by the points given above, it is capable of running FRG simulations from e.g.Density Functional Theory (DFT) or tight-binding (TB) datasets.The latter can be generated from within the code, with appropriate skeletons given to make usage simple and efficient.Furthermore the users are enabled to employ their own, custom momentum space meshings and thus can in principle use the code for conventional N -patch FRG simulations -even in the multi-orbital scenario.
We facilitate this by inclusion of an appropriate code skeleton.
The main computational challenge faced by this FRG implementation is memory size and bandwidth.The message passing interface (MPI) allows the vertex to be distributed on a large number of cluster nodes (with the restriction that splitting constrained to the bosonic momentum index q D ).This code is memory-bound; simulating e.g. a six band model on a 24 × 24 coarse momentum mesh would take 288 compute nodes requiring ∼ 140 GiB memory on each node to store the vertex objects.For fewer band models these requirements are drastically lowered and allow for quick and robust operation.

Code #2: TUFRG
The aim of the TUFRG implementation is to allow fast creation and testing of new material representations.This is enabled by removing unnecessary constraints imposed in the creation of previous frameworks, generalizing the implementation as much as possible, while extending to include SOC and multi-orbital systems.The cost introduced by this inclusion has necessitated a focus on performance, while the decision to generalize for all systems has necessitated some performance hits.
Due to the generalized nature of this implementation it will be worse in large unit cells compared to Code #3 and will be slower than Code #1 for grid-FRG.It is however an extremely accessible version of FRG and should be both easily usable and adaptable to uncharted problems.Furthermore it struggles much less under the memory constraining the calculations in Code #1, a similar calculation as described above would need a single compute node, MPI is required only for calculation speed.
The framework needs as input a tight-binding model with dispersion and interaction defined in momentum space as well as a momentum basis.From this the entirety of the formfactor expansion and TUFRG calculation will be performed.As the truncated unities are based heavily on grid-FRG the framework also provides a basic implementation of grid-FRG.While this frame-work is capable of doing the here published grid-FRG simulations, its primary aspiration is the TUFRG.

Code #3: RS-TUFRG
The objective of this implementation is to provide a flexible and easy to use real-and mixed space TUFRG algorithm.It is optimized for large unit cells, i.e. between 10 and a few 1000 sites per unit cell which allows for the study of many interesting phenomena, such as edge properties, disorder effects, behavior in quasicrystalline models and effects of different boundary conditions.The models also do not have to be defined on a lattice so that structures such as Barabasi-Albért networks can be analyzed.This broad application range of course leads to slight performance deterioration, thus it is advisable to resort to Code #1 or Code #2, for small unit cells or few band problems.By including all formfactor-bonds, this code can effectively perform the grid-RG calculations presented here.
The code itself is designed to be easily extendable with user defined models.For this it requires at least a definition of the underlying lattice or graph, a tight binding Hamiltonian, and a distance measure.It also allows for the inclusion of a single Matsubara frequency per channel during the flow.This enables the scaling test, where we check for the correct error scaling behavior of the interaction when compared to exact diagonalization, as shown in Appendix H.

Reproducibility: The Hubbard Model
To bolster confidence in the correctness of our results we first reproduce some established results of functional renormalization group calculations.The most basic model -thoroughly analyzed in Refs.[18, 19, 25, 28-30, 34-36, 39, 44-49] -is the Hubbard model for cuprates on a two-dimensional square lattice.We reproduce the results from Refs.[30,39], the next-nearest-neighbor hopping phase scan of the repulsive Hubbard model at the van Hove singularity.Note that the critical scale Λ c , which is derived from the artificial scale parameter Λ, is associated with the temperature of the phase transition and thus corresponds to a physical observable.Throughout this section, we therefore effectively compare both the critical temperature and the type of ordering to previous literature results.
The non-interacting part of the t−t Hubbard Hamiltonian is given by: The plots show the critical scale as a function of next-nearestneighbor hopping t with the chemical potential fixed at µ = −4t to pin the Fermi level to the van Hove singularity.Additionally, the type of instability is encoded as markers; circles: antiferromagnetism, plus signs: d-wave superconductivity, triangles: ferromagnetism.The interaction strength is set to U = 3.(a) Match of the non-TU results of Codes #1 and #2 to Ref. [30].(b) TUFRG results from Codes #2 and #3 matched to the results presented in Ref. [39].It should be noted that due to differing details of approximations and simulations we do not expect the results to exactly coincide.The discrepancies are pronounced in the proximity of the transition from superconducting to ferromagnetic phase at t ≈ −0.35.This region is both very dependent on implementational details as well as hard to properly resolve due to the low scales required (see Appendix B for the minimal scale allowed by our momentum resolution).This was already noted in the referred publications but is compounded by the different grid discretization chosen in (a) (grid-FRG vs N-Patch) leading to the behavior shown.
where the hopping amplitudes are t i,j = 1 if i, j are nearest neighbors and t i,j = t where i, j are nextnearest neighbors.We fix the chemical potential to µ = −4t pinning the system to van Hove filling.The interacting part of the Hamiltonian is governed by where n i,σ = c † i,σ c i,σ .Figure 3 demonstrates how all three codes presented in this work reproduce literature in terms of both the critical scales and the type of instability (details on the analysis of the effective vertex at the critical scale are presented in Appendix C).Note that our numerical implementations differ from each other and the reference in terms of the specifics and methods chosen within the scope of grid-FRG and TUFRG.To name examples, the type of integrator, cutoff scheme, stopping scale and termination criterion were not explicitly coordinated.We observe that the value of Λ c and the type of instability close to a phase transition are sensitive to these minor details, but overall the data shows agreement.The unattainable equivalence in Fig. 3 is a motivation for the exact comparison given in Section 5.The plots show the critical scale as a function of filling factor.Markers encode the type of ordering; circles: antiferromagnetism, plus signs: d-wave superconductivity.(a) Rashba model.The complete analysis can be found in Ref. [50], here we only show internal consistency between the three implementations.We use the following parameters: t = 0.15, α = 0.1, U = 3 and the filling factor ν ∈ [0.35, 0.59] with ν = 0 (ν = 1) corresponding to a completely empty (filled) system.(b) Graphene model.The markers and colors follow the same scheme as in (a).We choose the same parameter set as in Ref. [9]: t = 0.1, U = 3.6.For ease of comparison we use the doping δ = 2ν − 1 instead of the filling factor.The transition from magnetic to superconducting phase occurs at slightly different temperatures than in the reference which however employs singular mode FRG.Graphene is also notoriously resolution dependent (momenta and formfactors) in FRG calculations due to incommensurate phases [3] and the multi-site nature of the system.We nonetheless obtain the central results of the publication at δ ≈ 1/4 as well as the transition on either side.

Non-SU (2) Systems: The Rashba Model
To confirm all three codes' capabilities extend beyond the SU (2) symmetric one-band Hubbard model, we introduce a Rashba-z spin-orbit coupling (SOC) term to the kinetic part of the Hamiltonian breaking the SU (2) symmetry.It then reads where α is the SOC strength, σ is the vector of Pauli matrices, ij denote nearest neighbors and b ij are the corresponding directed bonds.Extensive coverage of Rashbaz SOC in the square lattice Hubbard model will be provided by a publication in preparation [50] where we employ FRG and the weak coupling renormalization group [51,52].Here we indicate only the results of a filling phase scan at fixed α = 0.1, t = −0.15 and U = 3 in Fig. 4 (a) to demonstrate internal consistency.

Multi-band Systems: Graphene
The third class of models we can check consistent results for are multi-site unit cell models.Ref. [9] provides singular mode FRG results for a tight-binding Hamiltonian of graphene to which we compare the three codes.The non-interacting part defined on the honeycomb lattice reads where o, o are the sites within a unit cell, i, j are unit cell indices and t ij,oo the hopping parameters.We set t ij,oo = 1 for nearest-neighbors and t ij,oo = t for nextnearest-neighbors.For the reproduction of the results in Fig. 4 (b) we choose t = 0.1.The site-dependent Hubbard interaction with U = 3.6 reads Additionally to the results in Fig. 4 (b), we have reproduced the critical interaction strength of the half-filled system (ν = 0.5) without next-nearest-neighbor hopping (t = 0) to be U crit ≈ 2.7t.This is in agreement with both Ref. [53] as well as the expectation that FRG lies between mean-field U crit ≈ 2.2t [54] and quantum Monte-Carlo U crit = 3.6t − 4.0t [55][56][57][58] predictions.

Results -Benchmark Systems
The test systems were chosen to cover a breadth of different aspects of functional renormalization group calculations in an attempt to reveal common errors.The two dimensional square lattice Hubbard model is chosen as the starting point.We choose the SU (2) symmetric representation which warrants the usage of the SU (2) symmetric set of flow equations from Fig. 1.To also cover the non-SU (2) flow equations it is imperative to include a non-SU (2) symmetric model.We therefore also consider the square-lattice Rashba model with spinorbit coupling.While this pair of models should complete the verification of the momentum dependencies in the contractions calculated during the flow, we want to extend the benchmarks to include a multi-band model with non-orthogonal lattice vectors.This is crucial in finding phase-errors as well as non-periodicities of the Hamiltonian.We therefore evaluate and publish data for a honeycomb lattice Hubbard model, which is a simple model for graphene.While this is obviously not a complete list of possible models one could show equivalence for, inconsistencies in this subset should uncover most inconsistencies that can arise during the implementation of the FRG.If more models are desirable or advisable the equivalence class will be extended and the repository [59] updated appropriately.

Square lattice SU (2): The Hubbard Model
The simplest test case is the square-lattice Hubbard model.Having indicated that our codes reproduce previously published results we now define parameters sets for the benchmarks.The chosen way of integration is grid-FRG -this allows the greatest confirmation across the three codes -other reproduction of data are available for comparisons upon request.Be mindful that RS-TUFRG or TUFRG benchmarks can only be supported by a subset of implementations.
Model: To make the test as accessible as possible we shall fix t = 1, this allows implementations which use this -rather common -scale to implement a comparative calculation.The remaining parameters of Eqs. ( 20) and ( 21) are (arbitrarily) defined to the following nonzero values: U = 3.0 and t = −0.1.In the case of the single band SU (2) symmetric Hubbard model the interaction is constant: The chemical potential is defined as µ = 0.5.We chose to define the chemical potential rather than the filling to remove the indirection of its calculation.This process is inaccurate, especially at the small system sizes discussed here.As the first step of any momentum space FRG calculation likely is checking for correctness in the noninteracting part of the Hamiltonian, we provide band structures along the high-symmetry paths for each of the models discussed in Appendix A.
Meshing: It is essential for exact numerical accuracy that the momentum space is meshed identically.We propose the following patching scheme: An equispaced grid of 6 points along each edge is laid within the first PZ to include all high-symmetry-points. 36 is a sensible choice for the total number of points because while being small enough to be quickly calculated it still contains points that are not high-symmetry points.The resulting mesh can be seen in Fig. 5 (a).This mesh is used wherever momenta are needed, for the discretization of the vertex as well as the integration over the loop momentum.It should be noted that the refinement mentioned in Section 2.4 is not applied in these calculations to reduce sources of potential error.
Flow: The most difficult aspect of the calculations to align are the parameters and calculations involving the integration from high to low scale.There are many options of integrators, Euler or adaptive, each with a significantly different update scheme for d/dΛ.Even minute details such as the order of operations (calculating ∆ Λ (Λ n+1 ) or ∆ Λ (Λ n )) are relevant.For this reason we decide to use In the hexagonal case we can see that the PZ mesh will fill the BZ only after backfolding of points, this is intended behavior.
the simplest possible integrator, constraining the issues in the numerical analysis: The integration scheme is a fixed-stepwidth Euler integration which yields We fix ∆ Λ = 0.1 ≡ const.and perform 90 iterations of the FRG-flow starting at Λ = 10.We disable all other termination conditions and are thus certain to obtain V (Λ = 1.0).To avoid misinterpretation, we provide pseudocode for the procedure: Algorithm 1 Flow scheme with constant update 1: Λ = 10 2: ∆ = 0.1 3: Initialize vertex V = V 0 4: for n = 0, 1, . . ., 90 do 5: Evaluate L(Λ) with Eq. ( 7) 6: Evaluate dV dΛ (Λ) via Fig. 1  7: V ← V + ∆ dV dΛ 8: For regulator of the FRG equations we also default to the simplest possible choice, the sharp cutoff.This is used in all subsequent calculations.

Square lattice non-SU (2): The Rashba Model
The meshing of the square BZ [Fig.5 (a)] as well as the parametrization of the Λ integration are equivalent to the setup presented for the Hubbard model.

Model:
We choose the model parameters to ensure broken particle-hole symmetry as well as broken SU (2) symmetry.The non-interacting Hamiltonian is defined by Eq. ( 22) while the interacting part Eq. ( 21) can be represented as a Hubbard interaction of electrons of opposite spins: We set the parameters to t = −0.1,U = 3.0, α = 0.5 and µ = 0.2.

Hexagonal lattice SU (2): The graphene Model
For the graphene model we choose the parametrization of the FRG-flow to remain equivalent but need to redefine the BZ meshing [Fig.5 (b)].

Model:
The parameters for the calculation of graphene are: U = 3, t = 0.1 and µ = 0.2.Due to the sublattice structure of graphene the choice of origin in the Fourier transforms of the Hamiltonian is relevant.This is detailed in Appendix D, here it suffices to mention that we choose the proper gauge which transforms both sites of the graphene model with respect to the origin of the unit cell.Any other choice of Fourier transform would result in an improper (non-periodic) gauge of the Hamiltonian.We specify the interaction in the orbital space of graphene to be a site-local Hubbard interaction: To be precise we specify the real space lattice vectors (a i ) and positions of the atoms within the first unit cell (x i ): Meshing: The meshing is defined as an equispaced grid of 6 × 6 points along the lines defined by the two reciprocal lattice vectors 0) T .This mesh covers the PZ [cf.Fig. 5 (b)] which is equivalent to any BZ due to the above choice of proper gauge.
Using this meshing ensures an even distribution of the weights to the momentum points during integrations.As before this mesh is used for all momenta required in the calculations.

TUFRG results
As Code #2 as well as Code #3 are capable of TUFRG calculations we are also able to generate benchmarks for this approximation.Equivalence is however much harder to achieve for TUFRG and because the target audience for the comparison is smaller the results are omitted from this publication.Upon request datasets as well as exact descriptions of the procedures involved will be provided.

Data Availability
The dataset is publicly available [59] for comparisons.Simulation codes will be made available from the corresponding authors upon reasonable request.They are expected to be made available in the near future.

Comparing Instructions
The comparative data includes the following aspects of the FRG evaluation: Discrepancies in these offer insight into possible faults in singular channels or prefactors.
The third and most important result to compare is the effective vertex at the final scale.Here we offer an array of values (of datatype double complex).The recommended approach is to sort the reference array as well as the respective result from the trial code by some metric (i.e.max[Re(z)] or max[Im(z)]) and compare the sorted arrays.Note that equivalence can be reached only to computational accuracy (discrepancies from the resolution of double, rounding, compiler optimizations, etc. may occur).
While it is also possible to obtain element-wise reproduction, this is dependent on the order of the discretized momentum points.Because this stems from the mesh generation scheme as well as the choice of PZ, neither of which has physical implications, we recommend subverting this dependency.

Conclusion and Outlook
This publication set out to rectify the missing link for internal consistency of momentum space functional renor-malization group calculations.We have verified the validity of the implementations by proving agreement with established results for momentum space functional renormalization group calculations and uniformity to unprecedented levels between the three implementations.This gives us the confidence in the claim that the datasets obtained are "correct" in the scope of the specific approximation (i.e.treating only the four-point vertex in momentum space while neglecting self-energies and frequency dependencies of the vertex) of FRG.We provide these results with the intent of their continued reproduction and verification by members of the FRG community.Because the realm of tests can be as vast as the scope of FRG calculations we invite a continued investment by the community.The obvious extensions include frequency and self-energy dependent calculations.Please engage with the authors for required assistance in the reproduction.
Furthermore, the authors are currently working to combine their codes under a single, versatile "community code" with a polished, common, easy-to-use interface.While this is a major undertaking it is necessary to make the resulting code accessible to both the FRG community and the more general audience of physicists interested in many-body phenomenae.
Fig. 6 Band structures of (a) square lattice tight binding model (Hubbard), (b) square lattice tight binding model with Rashba-SOC (Rashba) and (c) hexagonal lattice tight binding model (graphene), where the shorthand notation associates the band structures with the benchmark systems as defined in Section 5. We note for the Hubbard model the expected saddle point at the X-point.The square lattice Rashba band structure has band crossing points at X and M while for graphene we draw attention to Dirac cone around the K-point.

dl LΛ
+ (l, +(q − l)) (0, 0) Fig. 7 Integrated particle-hole (left) and particle-particle (right) loops dl LΛ ± (l, ±(q − l)) at high-symmetry values of q as a function of scale Λ.The points q are indicated in the legends.Note that at Λ = 10 −3 , we add a dashed gray emphasizing the minimum value of Λ for which the loops can be considered to be converged.Below this scale, the spacing of eigenenergies (for the discretization we chose) is larger than Λ leading to a divergence that is cut off.
For the superconducting phase we can consider the element at k 1 + k 2 = 0, using the remaining momentum dependencies to distinguish that we encounter a d-wave superconductor.In the truncated unities approximation this becomes even more powerful as the other momenta are already expanded into the formfactor basis.For the antiferromagnetic phase we consider the momentum transfer k 1 − k 3 = (π, π).As this drives the C-channel diagram we can assume a spin-densitywave with wave vector (π, π) which corresponds -on the square lattice -to antiferromagnetism.Similarly we can decipher ferromagnetism by examining the vertex at the momentum transfer k 1 − k 3 = (0, 0), a spin-densitywave with this wave vector corresponds to a ferromagnetic phase.

Susceptibilities
For particle-hole instabilities (spin-density waves and charge-density waves in the SU (2) case), it is instructive to study the (crossed) particle-hole susceptibility using the vertex given at the end of the FRG flow projected to the corresponding channels.[60] To efficiently calculate these, we define the Fermi particle-hole loop at scale Λ as with the Fermi function f (x) = 1 + e x −1 and X = b1 (k X ), X = b2 (k X ).From this, we calculate the fourpoint susceptibility as

Linearized Gap Equation
For P instabilities, it is convenient to solve a linearized gap equation to obtain details about the system's leading ordering tendencies.[60] We define the P -channel Fermi loop at scale Λ as We can then proceed to define a superconducting linearized gap equation as The eigenproblem in Eq. (C.4) is in general non-hermitian and thus numerically unstable.Therefore, we instead solve the following singular value decomposition: The right singular vectors V are gap functions projected to the Fermi surface (with "temperature" broadening set by Λ) and the left singular vectors do not include the Fermi surface projection and display the gap's symmetry.In the case of non-SU (2) and two-site systems, it is instructive to transform the superconducting (P ) gap to singlet [ψ(k)] and triplet [d(k)] space [61,62]: with σ Pauli matrices and σ0 the identity matrix.

Appendix D: Pitfalls during evaluations
During the production of the results for each model we encountered numerous possible pitfalls which can lead to different results in the test set.The first and most important one is to make sure that the models are defined exactly as above.Other common mistakes that occur are using different sign conventions or wrong input parameters, so if the results do not match this should be the first suspect.Otherwise the following section may help identify problems in the comparison.Most of these pertain to general problems which should be taken into account when implementing an FRG-code.

Signs and prefactors of the diagrams
One of the main issues with functional renormalization group calculations is fixing the prefactors (signs and values) within the flow equations (cf.Fig. 1).Here we want to list some strategies which can be employed in order to obtain the correct values for the special case of the square lattice Hubbard model at half filling (and t = 0).A graphical representation for the SU (2) symmetric rules is shown in Fig. 8.For all cases, we used the test system parameters, except for µ = 0, t = 0 such that the symmetries hold.

Non-SU (2) diagrams
The following rules apply to the three diagrams of the non-SU (2) equation: -For an initially attractive interaction, the absolute value maximum of the P -channel contribution must increase over the flow (for the first steps).
-Using an initially repulsive interaction, the C-and D-channel contributions to the effective vertex must increase over the flow (for the first steps).-The C-and D-channel contributions must be equivalent up to a reordering of the orbital and momentum indices and sign including the prefactors in the flow equations.

SU (2) diagrams
1.When calculating the single-channel flow in C-and P -diagrams the results must be equal when inverting U . 2. The total value of the D-channel contribution must be zero in the first step.3.For an initially repulsive interaction, the C-channel contribution to the effective vertex must increase over the flow (for the first steps).4. When calculating an attractive interaction in P -channel the maximum of the effective vertex must increase over the flow (for the first steps).

BZ-periodicity of systems with basis
Here we want to address the periodicity of the Hamiltonian in the case of existing sublattice structure.Fourier transforming into momentum space there are two options for handling the positions of the orbitals.We can either use the actual positions in the transformation or we map all sites onto the origin and Fourier transform with respect to this.Using the first option introduces non-periodicities.We will illustrate this using a simple two site model: We define the Fourier transforms with respect to the origin as follows: and resulting in the non-interacting Hamiltonian of the form which is a perfectly 2π periodic form of the non-interacting Hamiltonian.We shall call this gauge a "proper" gauge.Now to study the effects of an improper gauge, where we use the actual positioning of the sublattice within the Using these expressions to Fourier transform H 0 we can see that we collect some phase-terms which are not 2π periodic in k: We can therefore no longer evaluate H 0 after backfolding into the first BZ, instead we need to evaluate them in higher BZs.While in theory this is a seemingly natural gauge, the periodicity of the Hamiltonian is required for the momentum meshing used in this work.Furthermore using an improper gauge will change the topological properties of the model system.The definition of a 2π periodic gauged representation is therefore of paramount importance.

Handedness of Coordinate systems
While it may be obvious we want to briefly discuss the handedness of the momentum lattice.When defining a basis of three vectors spanning space we should define the third vector such that the span product v = (x×y)• z > 0. This is called a right-handed coordinate system.
When defining the real basis a i of the system instead we take care that it is right handed, the momentum lattice vectors G i can then be obtained via the transformations:

Appendix E: Refinement and Symmetries
Using the refinement to subsum the dependencies within the area associated with a mesh-point we run into issues which slightly break the symmetry of the Hamiltonian.
For clarity we will refer to the center the refinement collapses into as coarse point.

Definition of symmetric refinement
When defining the refined mesh points atop the coarse mesh points some attention is needed to obtain a set invariant under the symmetry operations of the coarse mesh.For an illustration of the set let us consider a hexagonal BZ into which we have introduced the coarse mesh via the methods discussed in Section 2.4.Defining the refined mesh as a scaled-down version of the coarse mesh would result in a mesh which is not invariant under the symmetry transformations, this is blatantly obvious from Fig. 10.The proper method employed in the implementations discussed here is to define a refined mesh of multiple times the size of the refined area size and reduce the points to the set which is closest to the coarse center point [Wigner-Seitz like construction, cf.Fig. 10 (b)].Depending on the system's lattice, it may be needed to double-count some of the refined points and introduce non-unitary integration weights.

Orbital-Band-Transformation symmetry breaking
Even when using a properly symmetric refined mesh there are remaining symmetry issues in the orbital-toband-transformations: We assume a model system with momentum-dependent orbital-orbital symmetry relations (such as the mapping of graphene-sites under rotation).
When we introduce refinement we have two options of calculating the averaged non-refined loop: -Calculate the average in band space, then transform into orbital space using the transformation matrices of the coarse points.This is subject to the bandcrossing problem presented in the paragraph "Band Crossing Problem" When refining in band space representation the band structure shown in the left plot we ignore the fact that a band crossing occurs within the region around the coarse mesh point.The orbital-to-band transformations, calculated only at the points of the coarse mesh point will map the averaged band into the orbital space respective their order at the coarse point.
-Transform each refined point into orbital space and average in orbital space.This is subject to the symmetry breaking of refinement presented the paragraph "Summed Symmetries".

Band Crossing Problem:
To illustrate the issue of evaluating the average in band space and transforming into orbital space afterwards we imagine a region of the BZ where multiple bands cross (cf.Fig. 11).Mathematically, we can formulate the problem that arises when the mapping from orbitals to bands changes as Summed Symmetries: Applying a symmetry to the Hamiltonian reads where U provides the needed phase-shifts and orbital maps for the transformation of orbitals across BZ-boundaries.
mixing of differing length scales in the interaction.This is exemplified in Fig. 12 where we can clearly see that the proper-gauge representation of the interaction (all sites Fourier-transformed with respect to the same position) allows the definition of the formfactor basis only with respect to the remaining A positions.The realspace representation of the first shell is drawn in Fig. 12 (a), where we also indicate an exemplary interaction from a B-site to all neighboring A-sites.If we now transform into the physical picture in (b) we can clearly see that the seemingly consistent interaction length mixes differing length scales.If we transformed into the formfactor basis the different length scales would have been integrated into the same formfactor shell, loosing the ability to differentiate their contributions.In the cross-channel projections this introduces a slight symmetry breaking into the approximation.We see breaking of otherwise degenerate states as well as slight modifications of critical scales.This behavior can be rectified by introducing a sitedependent formfactor expansion.We need to define a different set of neighboring formfactors for the A and B site in the above example.

Appendix G: Formfactor generation
We want to detail the generalization of the formfactor generation method described in Appendix 3 of Ref. [38].In order to do this we first recall the methodology of that publication before expanding to our new approach.

Original Method
Given a symmetry Group G of the problem with group elements g and representations Γ i with characters χ i (g) we define the projectors into the representation's contribution.
To now find the shell-basis for this representation apply the projector to a trial bond taken from the shell n of the realspace lattice φ n (r) = δ i,i+r .The result is the realspace-representation of the formfactor:  For the multi-dimensional representation we use a number of different trial bonds equal to the dimensionality.This is supposed to ensure that the number of formfactors is always sufficient to provide a unitary transformation from a realspace lattice shell into formfactor space.This however breaks down at longer ranges as will be seen in the next section.

Improved Version
The issues with the above mentioned method arise from its non-generality, the dimensionality of the representations is given but we encounter momentum shells which have more points than we have representations in the symmetry group (the fourth shell of the square lattice is an example).For these shells the above mentioned procedure does not generate a sufficient number of formfactors for the unitary transformation of the spaces.For these we are able to find multiple symmetry-inequivalent formfactors for a given representation.To produce these formfactors while maintaining orthogonality between the formfactors we use the following procedure: For all bonds of a given length in the realspace lattice we apply the projector defined above for all representations.We reduce the obtained (too large) set of realspace formfactors ϕ(r) to a linearly independent subset and orthogonalize using the Gram Schmidt orthogonalization procedure.
We can then promise the following equations for the formfactor basis: To ensure we generate the physical set of formfactors (the ones described in Ref. [38]) we use the following iterative algorithm for the generation of formfactors, this favors the formfactors generated by the original method.
Algorithm 2 Generate all formfactors of shell n The number of formfactors which need to be chosen is highly dependent upon the problem under scrutiny.If the interactions are highly localized, a few formfactors will suffice to capture the primary dependencies.A more detailed analysis of the necessary number of formfactors should be performed for each simulation.

Appendix H: Scaling tests
Beyond the equivalence checks and reproduction of known results, we further checked the correct scaling behavior of our implementations against exact solutions for small systems.Therefore a short exact diagonalization code for a 1D Hubbard chain, a 1D Rashba chain and 2D graphene has been implemented.We then calculated the occupation number ρ i,j = c † i c j using ED and a RS-TUFRG implementation in the single frequency approximation [27,41].The error is expected to scale proportionally to U 3 .To verify the correct scaling behavior, roughly 60 bosonic frequencies were necessary and the accepted error of the adaptive integrator has been set to 10 −6 , which allows us to verify the proportionality down to an absolute error of ≈ 10 −6 .
The results are summarized in Fig. 13.At low values of the interaction, the integration error leads to a deviation from the expected scaling, as can be seen in all three cases.For larger U, higher order terms can become important also deteriorating the scaling, this can be seen especially in the two 1D systems.Apart from these minor effects, the data follow the U 3 -curve very closely. Profe

Fig. 1
Fig. 1 Diagrammatic representation of the functional renormalization group equations.The lower set of diagrams represent the SU (2)-symmetric flow equations, where the SU (2)-degree of freedom is kept constant along the lines of the vertex.The line through the loop represents the scale derivative d/dΛ.We have already indicated the three distinct channels of the FRG flow, differentiated by the bosonic transfer momentum: P -, C-, and D-channel.

1 G 2 Fig. 5
Fig. 5 Momentum space meshing in test cases.(a) Square lattice models (Hubbard and Rashba) have a square Brillouin zone (BZ) and a square primitive zone (PZ) (b) Triagonal lattice systems (e.g.graphene) have a hexagonal BZ and a PZ in rhombus shape.Both meshes are an equispaced grid along the basis vectors of the reciprocal lattices G 1 and G 2 .In the hexagonal case we can see that the PZ mesh will fill the BZ only after backfolding of points, this is intended behavior.

1 .
The maximum contribution to ∆ Λ dV dΛ (Λ) by each of the P -,C-, D-channel diagrams (or the sum of the D-channel diagrams for the SU (2) systems) 2. The maximum vertex element after each iteration of the flow 3. The final effective vertex V (Λ = 1.0) at each discretization point of the Brillouin zone (and spin/site index for Rashba/Graphene) The first and second set of data can be compared by calculating the corresponding values in the trial implementation.The definition of maximum used is the maximal absolute value of the complex numbers |z| = Re(z) 2 + Im(z) 2 .

Fig. 8
Fig.8Graphical representation of test system properties.For all cases, we used the test system parameters, except for µ = 0, t = 0 such that the symmetries hold.(a) Vertex maximum as a function of scale in single channel flow in P for U = −3 (blue) and C for U = 3.The two lines are the same (point 1 in the text).(b) Full repulsive flow, D channel contribution that is identically zero at first step.(cf.point 2) (c) Full repulsive flow, all channels and vertex.(cf.point 3) (d) Full attractive flow (U = −3) with all channels and vertex.(cf.point 4) (a) Vertex maximum as a function of scale in single channel flow in P for U = −3 (blue) and C for U = 3.The two lines are the same (point 1 in the text).(b) Full repulsive flow, D channel contribution that is identically zero at first step.(cf.point 2) (c) Full repulsive flow, all channels and vertex.(cf.point 3) (d) Full attractive flow (U = −3) with all channels and vertex.(cf.point 4) unit cell to define the improper gauged Fourier transforms:

1 G 2 Fig. 10
Fig. 10 Refinement under symmetries.(a) Intuitive but incorrect definition of the refined mesh atop the coarse mesh.As can clearly be seen under a rotation of C 3 the indicated coarse points transform into each other, a rotation of the refinement however renders it incommensurate.(b) Proper definition of the refinement area as the Wigner-Seitz cell around the coarse point.It is apparent that these points will respect all symmetry transformations of the coarse mesh..