Benchmark calculations of multiloop pseudofermion fRG

Abstract

The pseudofermion functional renormalization group (pffRG) is a computational method for determining zero-temperature phase diagrams of frustrated quantum magnets. In a recent methodological advance, the commonly employed Katanin truncation of the flow equations was extended to include multiloop corrections, thereby capturing additional contributions from the three-particle vertex (Thoenniss et al. https://arxiv.org/abs/2011.01268; Kiese et al. https://arxiv.org/abs/2011.01269). This development has also stimulated significant progress in the numerical implementation of pffRG, allowing one to track the evolution of pseudofermion vertices under the renormalization group flow with unprecedented accuracy. However, cutting-edge solvers differ in their integration algorithms, heuristics to discretize Matsubara frequency grids, and more. To lend confidence in the numerical robustness of state-of-the-art multiloop pffRG codes, we present and compare results produced with two independently developed and algorithmically distinct solvers for Heisenberg models on three-dimensional lattice geometries. Using the cubic lattice Heisenberg (anti)ferromagnet with nearest and next-nearest neighbor interactions as a generic benchmark model, we find the two codes to quantitatively agree, often up to several orders of magnitude in digital precision, both on the level of spin-spin correlation functions and renormalized fermionic vertices for varying loop orders. These benchmark calculations further substantiate the usage of multiloop pffRG solvers to tackle unconventional forms of quantum magnetism.

Graphic abstract

1 Introduction

A fascinating phenomenon in the study of frustrated quantum magnets is the interplay of unconventional forms of magnetic order and the possible emergence of quantum spin liquid states near zero temperature [3]. The successful description of such low-energy states of quantum spin systems has, however, remained challenging, especially in the presence of competing interactions, geometric frustration, and in higher spatial dimensions.

Since its inception more than a decade ago [4], the pseudofermion functional renormalization group (pffRG) has become a powerful and flexible approach to map out the zero-temperature phase diagrams of various quantum spin models, both in two [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] and three spatial dimensions [16, 21,22,23,24,25,26,27,28,29]. Although the problem obtained after representing the spin operators by complex fermions is treated approximately, one of the striking features of pffRG is its ability to track competing instabilities in different interaction channels, allowing one to discriminate putative spin-liquid phases from long-range ordered magnetic ground states. This ability can be traced back [30, 31] to the inclusion of leading-order 1/S and 1/N diagrams (the former promoting classical magnetic order, the latter quantum fluctuations), which are treated on equal footing in pffRG by means of the routinely employed Katanin truncation [32].

Recently, the multiloop truncation scheme of the infinite hierarchy of fRG flow equations [33,34,35], previously used in the context of the Hubbard [36, 37] and Anderson impurity model [38], was applied to the zero-temperature pffRG by some of us [1, 2]. The convergence in the number of loops over a wide range of energy scales attested to the inner consistency of the pffRG method, despite being used in the strong-coupling limit. These developments were accompanied and facilitated by substantial improvements of the numerical implementation that remedy many shortcomings of previous studies. Yet, some of these advances, such as the employed integration routines and adaptive Matsubara frequency grids [1, 2], rely on certain numerical heuristics, affecting, e.g., the minimal grid spacing and largest Matsubara frequencies considered. Therefore, quantitative agreement between different implementations is, although highly desired, not guaranteed a priori.

In the present work, we provide evidence for the numerical robustness of pffRG by benchmarking two independent state-of-the-art solvers, one provided by a research group at LMU Munich (dubbed code #1 in the following), and one by a Cologne–Würzburg collaboration (denoted by code #2) with an open-source release [39]. As a test case, we consider ferro- and antiferromagnetic Heisenberg models on the simple cubic lattice and compare our results both on the level of renormalized couplings (i.e. fermionic vertex functions) as well as for the (post-processed) spin-spin correlation functions.

The remainder of the paper is structured as follows. We begin by providing a brief overview of the multiloop pffRG in Sect. 2. This is followed by an in-depth comparison of the numerical results produced by the two codes at hand in Sect. 3. Finally, in Sect. 4, technical aspects of the implementation, such as the choice of frequency grids, integration routines and differential equation solvers are discussed, with special emphasis devoted to their influence on the numerical stability and accuracy of the two codes.

2 Multiloop pseudofermion fRG

Within the pffRG approach, one can study generic spin-1/2 Hamiltonians with bilinear spin couplings, i.e.,

$$\begin{aligned} \mathcal {H} = \tfrac{1}{2} \sum _{ij} J^{\mu \nu }_{ij} S^{\mu }_i S^{\nu }_j \,. \end{aligned}$$

(1)

Here, the spin operators $S^{\mu }_i$ live on the sites i of an arbitrary lattice, and the exchange matrices $J^{\mu \nu }_{ij}$ are assumed to be real. The spin operators are represented in terms of complex pseudofermions $f^{(\dagger )}_{i \alpha }$ with $\alpha \in \{\uparrow , \downarrow \}$ as

$$\begin{aligned} S^{\mu }_i = \tfrac{1}{2} \sum _{\alpha , \beta } f^{\dagger }_{i \alpha } \sigma ^{\mu }_{\alpha \beta } f_{i \beta } \,, \end{aligned}$$

(2)

where $\sigma ^{\mu }_{\alpha \beta }$ for $\mu \in \{x, y, z\}$ are the Pauli matrices. This yields a purely quartic Hamiltonian which can be treated by established functional RG techniques.

Note that the pseudofermion representation of the spin algebra comes with an artificial enlargement of the local Hilbert space dimension, which must be dealt with by an additional particle number constraint $\sum _{\alpha } f^{\dagger }_{i\alpha } f_{i \alpha } = 1$ on every lattice site. In practice, this constraint is not enforced, but holds on average due to particle-hole symmetry [1, 2, 4]. Nevertheless, the influence of fluctuations can be quantitatively gauged by explicitly computing the variance of the number operator, which can be expressed through the equal-time spin-spin correlation function $\langle S_{i}^{\mu } S_{i}^{\mu } \rangle$ [1]. Although fluctuations are not fully suppressed, even if a local level repulsion term $A S^{\mu }_i S^{\mu }_i$ (with $A < 0$) is employed, recent studies [1, 19, 23, 30] pointed out that observables extracted from pffRG flows are qualitatively unaffected by the unphysical Hilbert space sectors.

An alternate decomposition of the spin operators into Majorana instead of Abrikosov fermions allows one to circumvent the problem of unphysical states in the fermionic representation at the cost of redundant copies of physical Hilbert-space sectors [40]. For moderately high temperatures, the latter approach was recently shown to enable an accurate calculation of thermodynamic observables [41], such as the free energy and specific heat. However, the approach was also found to suffer from unphysical divergencies when approaching the $T \rightarrow 0$ limit, which we consider here (for the Abrikosov fermion decomposition).

Since kinetic contributions are absent in the pseudo-sfermion representation of Eq. (1), the free propagator assumes the simple form

$$\begin{aligned} G_0(1'|1)&= (i \omega _{1})^{-1} \delta _{i_{1'} i_{1}} \delta _{\alpha _{1'} \alpha _{1}} \delta (\omega _{1'} - \omega _{1}) \,, \end{aligned}$$

(3)

diagonal in all indices. To successively integrate out high-energy modes and thus provide an effective low-energy description of a given model, a cutoff parameter, here denoted as $\Lambda$, is introduced in the bare propagator. The fRG equations then govern the flow of the n-particle vertices from the UV limit $\Lambda \rightarrow \infty$, where the regularized bare propagator vanishes, to the infrared limit $\Lambda \rightarrow 0$, where one recovers the physical theory. As such, there is a certain degree of freedom in the cutoff implementation. A popular choice for the regulator in pffRG is a Heavyside step function, which sharply suppresses frequency contributions $|\omega | < \Lambda$. This choice is very useful for analytical treatments of pffRG in the large-S and large-N limit, where the flow equations can be solved exactly and reproduce mean-field gap equations [30, 31]. However, if numerical calculations are employed away from these limits, a non-analytic regulator spoils the smoothness of the right-hand side of the flow equations, and therefore limits the applicability of higher-order integration routines. For this reason, we consider a smooth regulator

$$\begin{aligned} R^{\Lambda }(\omega ) = 1 - e^{-\omega ^2 / \Lambda ^2} \,, \end{aligned}$$

(4)

throughout this manuscript, and implement the cutoff as $G^{\Lambda }_0(\omega ) = R^{\Lambda }(\omega ) G_0(\omega )$, with $G_0(\omega ) \equiv (i \omega )^{-1}$.

To make the infinite hierarchy of fRG flow equations amenable to further calculations, a truncation is necessary. Usually, this is done by neglecting all n-particle vertices of $n=3$ and higher [32]. However, to capture the physics of interest in pffRG, one must already go beyond that using the Katanin truncation, which feeds the $\Lambda$ derivative of the self-energy $\Sigma ^{\Lambda }$ back into the flow of the two-particle vertex $\Gamma ^{\Lambda }$ [4]. Within this truncation, the flow equations schematically read

$$\begin{aligned} \frac{{d}}{d\Lambda }\Sigma ^\Lambda&= - \big [\Gamma ^{\Lambda } \circ S^{\Lambda } \big ]_{\Sigma } \,, \end{aligned}$$

(5)

$$\begin{aligned} \frac{{d}}{d\Lambda }\Gamma ^\Lambda&= \sum _c \dot{\gamma }_{c}^\Lambda = -\sum _c [\Gamma ^\Lambda \circ \partial _\Lambda (G^\Lambda \times G^\Lambda ) \circ \Gamma ^\Lambda ]_{c} \,. \end{aligned}$$

(6)

Here, we introduced the loop function $[\Gamma \circ G]_{\Sigma }$ and the single-scale propagator $S^{\Lambda } \equiv -\frac{{d}}{d\Lambda }G^{\Lambda }|_{\Sigma ^{\Lambda } = \text {const.}}$. We categorized the contributions to the flow of $\Gamma$ into three distinct channels c: the particle-particle (s) channel, the direct particle-hole (t) channel, and the crossed particle-hole (u) channel. Each “bubble” term, with the general form $[\Gamma \circ (G \times G') \circ \Gamma ']_c$, describes the flow of a two-particle reducible vertex $\gamma _c$. As all self-energies, vertices, and related correlators are $\Lambda$-dependent, we refrain from writing this dependence explicitly in the following.

The multiloop fRG (mfRG) flow [33,34,35], recently employed within pffRG [1, 2], is an attempt to go beyond the Katanin truncation and capture even more contributions from n-particle vertices with $n\ge 3$. It can be derived from the parquet approximation [42], which self-consistently connects one- and two-particle correlation functions via the Schwinger–Dyson (SDE) and Bethe–Salpeter equations (BSE), and as such the inherent dependence of the $\Lambda \rightarrow 0$ fRG result on the specific choice of regulator is eliminated [34]. This approximation includes all those contributions to the flow of the two-particle vertex which can be efficiently calculated, i.e., with the same cost as the one-loop flow in Eqs. (5) and (6). Summarized briefly: To obtain the mfRG flow of $\gamma _c$, one iteratively computes multiloop corrections to the one-loop ($\ell =1$) result, using bubble functions with undifferentiated propagators but differentiated vertices. In a similar fashion, one can recover equivalence to the SDE, by feeding back the so-determined vertex corrections into the self-energy flow.

One of the most important ingredients to achieve sufficient numerical accuracy throughout the multiloop flow is an appropriate treatment of the frequency dependence of the two-particle vertex. In Ref. [43], a parametrization in terms of one bosonic and two fermionic frequencies (the fourth frequency argument is fixed by energy conservation) for each two-particle reducible vertex was put forward. This parametrization captures the non-trivial high frequency asymptotics of the vertices while being numerically efficient. Code #1 uses precisely the proposal of Ref. [43], and the diagrams contributing to each channel are grouped into four asymptotic classes $K_n$ as

$$\begin{aligned} \gamma _{c}(\omega _c, \nu _c, \nu '_c)&= K_{1, c}(\omega _c) \nonumber \\&\quad + K_{2, c}(\omega _c, \nu _c) + K_{2', c}(\omega _c, \nu '_c) \nonumber \\&\quad + K_{3, c}(\omega _c, \nu _c, \nu '_c) \,, \end{aligned}$$

(7)

where we displayed only frequency arguments for brevity. Here, $\omega _c, \nu _c$ and $\nu '_{c}$, denote the natural frequency arguments for diagrams reducible in channel c (see Ref. [1] for the conventions used). The $K_{n}$ asymptotically decay to zero in each frequency, allowing one to reduce the necessary number of arguments when summing up the asymptotic classes to obtain $\gamma _c$. Code #2 chooses a slightly different approach, by defining asymptotic classes $Q_n$ [44] as

$$\begin{aligned} Q_{1, c}(\omega _c)&= K_{1, c}(\omega _c) \nonumber \\ Q_{2, c}(\omega _c, \nu _c)&= K_{1, c}(\omega _c) + K_{2, c}(\omega _c, \nu _c)\nonumber \\ Q_{2', c}(\omega _c, \nu '_c)&= K_{1, c}(\omega _c) + K_{2', c}(\omega _c, \nu '_c)\nonumber \\ Q_{3, c}(\omega _c, \nu _c, \nu '_c)&= K_{1, c}(\omega _c) \nonumber \\&\quad + K_{2, c}(\omega _c, \nu _c) + K_{2', c}(\omega _c, \nu '_c) \nonumber \\&\quad + K_{3, c}(\omega _c, \nu _c, \nu '_c) \,, \end{aligned}$$

(8)

with the respective choice of natural frequency arguments outlined in Ref. [2]. Since the $K_n$ classes decay to zero for large frequencies, the $Q_n$ (at least for $n > 1$) are projected to a lower class. For instance, $Q_{3, c}(\omega _c, \nu _c, \nu '_c) = Q_{2, c}(\omega _c, \nu _c)$ if $|\nu '_c| \rightarrow \infty$. Let us emphasize that both parametrizations contain the same information about the asymptotic structure of the two-particle vertices, as the $K_n$ and $Q_n$ parametrizations can be exactly transformed into each other. For an appropriate choice of numerical frequency grids, both parametrizations are therefore equally valid and differ only in numerical performance. The former approach allows for a more fine-grained adjustment of discrete frequencies to the asymptotic decay of individual classes, while the latter reduces the cost of evoking a two-particle vertex from a summation of up to four classes $K_n$ to loading just a single $Q_n$.

The central observable computed from the pffRG equations is the flowing spin-spin correlation function,

$$\begin{aligned} \chi ^{\mu \nu }_{ij}(i \omega = 0) = \int _{0}^{\infty } d\tau \langle T_{\tau } S^{\mu }_{i}(\tau ) S^{\nu }_{j}(0) \rangle \,, \end{aligned}$$

(9)

where we omit indication of the $\Lambda$-dependence for brevity. In all models considered here, the interactions in the Hamiltonian are diagonal and $\mathrm {SU}(2)$-symmetric. This leads to spin-spin correlations that are symmetric as well, and we thus define $\chi _{ij} \equiv \chi ^{xx}_{ij} = \chi ^{yy}_{ij} = \chi ^{zz}_{ij}$.

The spin-spin correlations can be used to identify transitions into phases with broken symmetries; there, the flow becomes unstable at some $\Lambda _{\mathrm {T}}$ and must be stopped. For long-range ordered states, the momentum ${\varvec{k}}$ for which the structure factor

$$\begin{aligned} \chi ({\varvec{k}}, i\omega ) = \frac{1}{N_{\text {sites}}} \sum _{ij} e^{i {\varvec{k}} \cdot ({\varvec{R}}_i - {\varvec{R}}_j)} \chi _{ij} (i\omega ) \end{aligned}$$

(10)

(i.e. the Fourier transform of $\chi _{ij}$) is most dominant gives an indication of the emergent magnetic order, as exemplified in Fig. 1. A smooth flow down to the infrared $\Lambda \rightarrow 0$ is, on the other hand, associated with non-magnetic phases, such as spin liquids, dimerized, or plaquette-ordered states.

3 Results

Fig. 1

Momentum-resolved structure factors within the first Brillouin zone of the cubic lattice for (a, b) the ferromagnetic case at $\Lambda /J = 0.8$ and (c, d) the paramagnetic case at $\Lambda /J = 0.3$, computed for (a, c) $\ell = 1$ and (b, d) $\ell = 3$ using code #2. The ferromagnet shows a sharp peak at the $\mathbf {\varGamma }$ point, without visible difference between the two loop orders. The putative paramagnet shows a broadened distribution of spectral weight centered around soft maxima at the ${{\varvec{M}}}$ points in $\ell =1$ calculations, while the structure factor peaks more distinctively for $\ell =3$, signalling the onset of magnetic order instead

Fig. 2

Inverse spin-spin correlation function for the ferromagnet as a function of $\Lambda$. Shown here is a comparison of the $\ell =1$ and $\ell =3$ flows obtained from both codes. The dotted line is a $\Lambda ^{-1}$ fit [$\chi _{\mathrm {C}}= CJ/(\Lambda - \Lambda _{\mathrm {C}})$] to the data at $\Lambda /J\in \left[ 1.0, 4.0\right]$. The transition to a ferromagnetically ordered phase is visible as a sharp downturn away from Curie–Weiss behavior. Inset: Definition of the first, second, and third nearest-neighbor interaction, $J_1$ (green), $J_2$ (purple), and $J_3$ (yellow)

To benchmark the two codes, we calculate the spin-spin correlations and pseudofermion vertices of an extended Heisenberg model on the cubic lattice with a maximum correlation length $\xi = 5$ in units of the lattice spacing [1]. The corresponding three-dimensional cluster contains $N = 515$ sites, small enough to efficiently compare the two codes but large enough to produce the (qualitatively) correct physics. The corresponding Hamiltonian with up to third-neighbor interactions (see inset in Fig. 2) reads

$$\begin{aligned} \mathcal {H} = J_1 \sum _{\langle ij\rangle } S^{\mu }_i S^{\mu }_j + J_2 \sum _{\langle \langle ij\rangle \rangle } S^{\mu }_i S^{\mu }_j + J_3 \sum _{\langle \langle \langle ij\rangle \rangle \rangle } S^{\mu }_i S^{\mu }_j \,,\nonumber \\ \end{aligned}$$

(11)

where we fix $J \equiv \sqrt{J_1^2 + J_2^2 + J_3^2}$ as the unit of energy. We focus on two choices of these interaction parameters to highlight differences between fRG flows in different phases:

$$\begin{aligned}&J_1< 0, \qquad \qquad J_2 = 0, \qquad \qquad J_3 = 0, \end{aligned}$$

(12)

$$\begin{aligned}&J_1> 0, \qquad \quad \!\!\! J_2/J_1= 0.6, \qquad \! J_3/J_1= 0.25, \end{aligned}$$

(13)

where Eq. (12) yields a nearest-neighbor ferromagnet and the setup of Eq. (13) was previously reported to result in a paramagnetic ground state [21].

Rewriting each spin operator $S^\mu$ in the Hamiltonian in terms of pseudofermions leads to an expression proportional to ${f^\dagger _{\alpha '}} {f_{\alpha }} {f^\dagger _{\beta '}} {f_{\beta }}$, with interactions proportional to $\sum _\mu \sigma ^{\mu }_{\alpha '\alpha } \sigma ^{\mu }_{\beta '\beta }$. Exploiting this $\mathrm {SU}(2)$ symmetry (the interactions are diagonal and of equal magnitude in every spin direction), the flowing pseudofermion vertex $\Gamma$ (and each of its two-particle reducible parts $\gamma _c$) can be decomposed into a spin component $\Gamma ^s$, proportional to the latter combination of Pauli matrices, and a density component $\Gamma ^d$ proportional to $\delta _{\alpha '\alpha }\delta _{\beta '\beta }$ [4, 45]. Note that the density component, although initially vanishing for any typical spin model, becomes finite away from the UV limit and is essential for tracking the evolution of all symmetry-allowed couplings under the RG flow.

Fig. 3

Structure factor for the ferromagnet along a high-symmetry path of the cubic lattice Brillouin zone. The results are in excellent agreement between both codes, both for $\ell = 1$ and $\ell = 3$, showing dominant ferromagnetic correlations indicated by a sharp peak around the $\mathbf {\varGamma }$ point. Inset: Zoom into the path segment connecting the ${\varvec{X}}, {\varvec{M}}$, and ${\varvec{R}}$ point

Fig. 4

Frequency structure of self-energy and $t$-reducible vertex for the ferromagnet at different values of $\Lambda / J$ for $\ell = 3$ flows. The self-energy is purely imaginary and antisymmetric in frequency space, while all vertex components are real and symmetric along the directions plotted here. We show two cuts through the three-dimensional structure of $\smash {\gamma _{t, \langle ij\rangle }^{\Lambda , \mu } (\omega , \nu , \nu ')}$: A cut along the bosonic frequency axis $\omega$, with both fermionic frequencies set to $\nu = \nu ' = 0$, and a cut with equal fermionic frequencies $\nu = \nu '$, where the bosonic frequency was set to $\omega = 0$. The first cut is not shown for $\gamma _t^d$ as $\smash {\gamma _{t, \langle ij\rangle }^{d}(\omega , 0, 0) = 0}$ due to symmetry [1, 2]. The most prominent structure in the $t$-reducible vertex is a peak around zero bosonic frequency $\omega = 0$ that grows in magnitude and becomes sharper as $\Lambda$ is decreased. This indicates ferromagnetic correlations that grow stronger as the ordering phase transition is approached. In all components, there is quantitative agreement between the two codes

3.1 Ferromagnetic phase

With pure nearest-neighbor ferromagnetic interactions, the zero-temperature ground state is intuitively expected to be a ferromagnet. Therefore, in the context of pseudofermion fRG, there should be a transition at some finite $\Lambda _{\mathrm {T}}> 0$ from a paramagnetic regime at large $\Lambda > \Lambda _{\mathrm {T}}$ to the ferromagnetic phase at $\Lambda < \Lambda _{\mathrm {T}}$. Approaching the transition, the spin-spin correlator $\chi _{ij}$ is expected to diverge, similar to a finite-temperature phase transition. In this case, a peak will form at the $\mathbf {\varGamma }$ point in reciprocal space, as is visible Fig. 1, since the correlations are uniform and positive in a ferromagnet.

Close to the transition, the flow is supposed to visibly deviate from its paramagnetic Curie–Weiss behavior $\chi _{ii} \approx {CJ}/{(\Lambda - \Lambda _{\mathrm {C}})}$ at large $\Lambda \gg \Lambda _{\mathrm {T}}$. For this reason, it is convenient to plot the inverse correlator $1/\chi _{ii}$ as a function of $\Lambda$ to locate the transition, as shown in Fig. 2. Here, the $1 / \Lambda$ behavior appears as a straight line with slope $1/C$ displaced horizontally by $\Lambda _{\mathrm {C}}/J$ and the transition to the ferromagnetic phase is visible as a sharp turn down to a smaller inverse correlation function at $\Lambda /J \approx 0.76$. The structure factor at $\Lambda$ close to $\Lambda _{\mathrm {T}}$, shown in Figs. 1 and 3, has a single peak at the $\mathbf {\varGamma }$ point, signifying an instability towards ferromagnetic order. This, as well as the Curie–Weiss fit parameters, are consistent across both considered loop orders $\ell = 1, 3$ and both codes, while $\Lambda _{\mathrm {T}}$ differs slightly.

Since both implementations obtain the spin-spin correlations by post-processing the vertices, any discrepancy therein originates from differences in the vertices. Therefore, a more detailed examination of the $1 / \chi _{ii}$-deviations between the codes for $\ell = 1$ will follow once the flow of the vertex components has been discussed. Moreover, even if the flows for the $\chi _{ij}$ agree perfectly (as, e.g., in the regime $\Lambda > \Lambda _{\mathrm {T}}$), discrepancies in the vertices cannot be fully excluded, as post-processing spin-spin correlations from pseudofermion vertex data amounts to integrating a combination of several propagators and the vertex over two frequencies [1]. Hence, this additional step might hide potential differences in the vertex data.

Fig. 5

Decomposition of the $\smash {\gamma _{t, \langle ij\rangle }^{s}(\omega , \nu , \nu ')}$ vertex for the ferromagnet into asymptotic classes $K_{1,t}, K_{2,t}, K_{3,t}$ (first, second, third row) for the $\ell \!=\! 3$ flows at $\Lambda / J \!=\! 0.8$. Frequency axes shown here are the same as in Fig. 4. As the flow is close to the ordering phase transition at this value of $\Lambda$, strong ferromagnetic correlations are present as a peak around $\omega = 0$ in $K_{1,t}$. The other classes are at least one order of magnitude smaller. In all classes, both codes show quantitative agreement

To investigate this further, we focus on the $t$-reducible vertex $\gamma _t$ plotted in Fig. 4 at various values of $\Lambda$: Its spin component $\gamma _t^{s}$ (second and third column) is responsible for the transition and becomes sharply peaked at small bosonic frequencies $\omega \approx 0$. Its density component $\gamma _t^{d}$ (last column) with its extended structures and peaks at non-zero fermionic frequencies $\nu$ is particularly difficult to resolve and thus most likely to contain numerical artifacts. Comparing $\gamma _t$, as well as the the self-energy $\Sigma$ between the codes, we find quantitative agreement also on this very detailed level of inspection.

As outlined in Sect. 2, both codes use a decomposition of the reducible vertices $\gamma _s, \gamma _t, \gamma _u$ into four asymptotic classes each. The decomposition into asymptotic classes $K_n$ is shown for $\gamma _t^{s}$ at $\Lambda / J =0.8$ in Fig. 5, where we omit $K^s_{2',t}$, as it is equal to $K^s_{2,t}$ by crossing symmetry [1, 2]. Note that, while these vertices can directly be extracted from code #1, an additional transformation is applied to the $Q_n$ decomposition of code #2 [see Eq. (8)]. The peak in $\gamma ^{s}_{t}$ at small bosonic frequencies in Fig. 4 is found to stem from the $K_1$ contribution, which is an order of magnitude larger than the other classes. In $K_2$ and $K_3$, extended structures with multiple maxima and minima exist. It is thus crucial to use a frequency mesh with enough mesh points in an extended region around the origin to control numerical interpolation errors (see Sect. 4).

Though the codes implement the vertex decomposition differently (see Sect. 2) and use different approaches to build appropriate frequency meshes (see [1, 2] for a detailed description), all components of the vertex are consistent with each other. This demonstrates that it is possible to gain control over said interpolation errors by a careful adaptive implementation that places enough mesh points where they are needed.

Since the numerical error incurred by interpolation of the continuous frequency structure from a discrete mesh is particularly relevant whenever sharp structures are present in the vertex, different choices of frequency meshes have strong effects close to phase transitions, where some couplings are expected to diverge. For instance, in the ferromagnetic setup discussed above, the transition was induced by a peak in the spin component of the $t$-reducible vertex that grows quickly and starts to diverge, as can be seen in the second column of Fig. 4. As the transition is approached, this peak progressively becomes sharper and thus more difficult to resolve using discrete meshes. Thus, minor differences in mesh spacing can induce differences in the flow at the transition, though the qualitative, physical results remain unchanged.

To investigate the effects of changes in the mesh spacing explicitly, we compared results obtained from both codes with artificially modified meshes. Both implementations make use of adaptive frequency grids where, during the flow, the mesh spacing is adjusted according to the frequency structure of the vertex. The simplest way to manipulate the meshes is to rescale them by an artificial scaling factor $\kappa$. In Fig. 6, we show the effect of such a rescaling on the $\ell =1$ flow from Fig. 2. Above $\Lambda /J \approx 0.8$, all frequency structures in the vertex are fairly broad and easy to resolve. Consequently, rescaling the frequency grid has little effect and values $\kappa = 0.5 \ldots 3.0$ result in the same flow and also the same Curie–Weiss fit parameters. Below that point, the flows differ more and more as structures become sharper and ultimately predict slightly different transition points $\Lambda _{\mathrm {T}}/J$. Nevertheless, all flows predict a transition to the same ferromagnetic phase, which can be identified by a peak in the structure factor at the $\mathbf {\varGamma }$ point.

Fig. 6

Flows with rescaled frequency meshes. Comparison of the flow of inverse static on-site spin correlations $1/\chi _{ii}(i\omega =0)$ obtained using frequency meshes with different scaling factors $\kappa$. The dotted line is a $\Lambda ^{-1}$ fit to the data at $\Lambda / J \in [1.0, 4.0]$. For all values of $\kappa$, a transition to a ferromagnet is visible as a sharp turn down. The predicted transition point as well as the slope of $\chi$ in the region $\Lambda /J < 0.8$ differs, while the behavior at large $\Lambda > J$ remains identical

Fig. 7

Inverse spin-spin correlation function for the putative paramagnet as a function of $\Lambda$. Shown here is a comparison of the $\ell =1$ and $\ell =3$ flow obtained from both codes. The dotted line is a fit of a $\Lambda ^{-1}$ power law to the data at $\Lambda / J \in [1.0, 4.0]$. For $\Lambda /J \ge 0.5$, the $\Lambda ^{-1}$ behavior is followed almost perfectly. At smaller $\Lambda /J$, the $\ell =1$ and $\ell =3$ flows disagree: The $\ell =1$ curve smoothly approaches $\Lambda =0$ (staying above the power law), indicating antiferromagnetic correlations. By contrast, the $\ell =3$ curve displays a downward cusp, similar to Fig. 2, and thus predicts an ordered state

3.2 Paramagnetic phase

For the second set of parameters, Eq. (13), all interactions up to the third neighbor are antiferromagnetic. Consistent with prior work using one-loop fRG [21], both codes find a paramagnetic ground state for $\ell =1$, indicated by a smooth and regular flow down to $\Lambda =0$ in Fig. 7.

Remarkably, the $\ell =3$ data predicts a qualitatively different phase: There is a divergence in the spin correlations at $\Lambda _{\mathrm {T}}/ J \approx 0.24$, indicating an ordering transition at a scale roughly three times lower than for the ferromagnetic ordering instability discussed in the previous section. Such a reduced ordering scale is not unexpected for an exchange-frustrated spin system when compared to an unfrustrated one, but sometimes hard to establish.

Probing the structure factor in the vicinity of the divergence reveals a strong enhancement of magnetic correlations compared to the $\ell = 1$ flow, as indicated by sharpened Bragg peaks around the ${\varvec{M}} = (0, \pi , \pi )$ points in Figs. 1 and 9. These correspond to antiferromagnetic correlations between planes orthogonal to the vector connecting the second nearest-neighbors along diagonals of the faces in the cubic unit cell (shown in purple in Fig. 2). Our result is consistent with earlier observations of long-range $(0, \pi , \pi )$ order neighboring the paramagnetic phase [21]. Yet, the mfRG flows obtained from both codes suggest a rather strong modification of the respective phase boundaries as the coupling parameters investigated here were previously predicted to be deep in the non-magnetic regime.

Fig. 8

Frequency structure of self-energy and $t$-reducible vertex for the putative paramagnet at different values of $\Lambda / J$ for $\ell = 1$ and $3$ flows. As the $\ell = 3$ flow diverges at $\Lambda / J \approx 0.24$, only $\ell = 1$ is shown at $\Lambda / J = 0.05$. The same cuts through the three-dimensional frequency structure of the vertices are shown as in Fig. 4. Again, a peak in the $\gamma ^{s}_{t, \langle ij\rangle }$ component (second column) indicates strong correlations that become stronger as $\Lambda$ is further decreased. In contrast to the ferromagnetic case, this peak is negative, indicative of antiferromagnetic correlations, and there is a sizeable contribution of $\gamma ^s_t$ for nonzero fermionic frequencies $\nu , \nu '$ (third column), particularly for $\ell = 3$

In the vertex (see Fig. 8) and self-energy, there is again very good quantitative agreement between both codes. At $\Lambda /J = 0.05$, small quantitative differences between code #1 and #2 appear in the density component $\gamma _t^d$ of the $t$-reducible vertex, consistent with the earlier remark that it is the most difficult component to resolve well.

The $\ell = 1$ and $\ell = 3$ flows are very similar down to $\Lambda /J \ge 1$. Contributions of $\ell > 1$ terms become significant at $\Lambda /J \approx 1$ and eventually lead to an ordering instability induced by a peak in the $\smash {\gamma _{t}^{s}}$ component that diverges at $\Lambda / J \approx 0.24$. In contrast to the ferromagnetic case, this peak is negative, indicating anti-correlation. Along the fermionic $\nu$ frequency axis, the vertex shows an extended structure with multiple peaks of similar magnitude to the one on the bosonic axis. Since the $K_1$ class has no fermionic frequency, this means that, remarkably, other classes reach an order of magnitude comparable to $K_1$, as shown explicitly in Fig. 10. Consequently, vertex structures along fermionic frequency axes, in contrast to the ferromagnetic transition, become sizeable. It is therefore crucial to resolve the full three-dimensional frequency structure in $K_3$. Though numerically expensive, a large number of mesh points is necessary to ensure sufficient accuracy, as inadequate resolution of features along the fermionic frequency axes can strongly affect the fRG flow. This is even more important for multiloop flows, where interpolation errors might accumulate during the iteration over loop orders.

4 Technical aspects

To conclude our benchmark calculations, we discuss some of the particularly relevant technical aspects (see Table 1) which are needed to obtain confidence that we have sufficient degree of control over numerical errors. In doing so, we will also connect to the existing literature and scrutinize some of the algorithmic approaches which are routinely employed in the pffRG community.

4.1 Frequency grids

Both the self-energy and two-particle vertices are functions of Matsubara frequencies, which are continuous in the zero-temperature limit. A numerical implementation has to sample these functions on a finite grid and interpolate their values in between the sampling points. In many previous works (see e.g. Refs. [4, 19, 46]), the same frequency grid was chosen for the self-energy and all reducible vertices, usually featuring logarithmically increasing distances between adjacent grid points starting from some small but finite frequency. The intention behind such a choice of frequencies was to resolve the structure around zero frequency with high accuracy while coarse-graining high-frequency tails. Moreover, each vertex component was parametrized in terms of the three bosonic transfer frequencies, instead of the channel-specific mixed bosonic-fermionic frequency treatment utilized by codes #1 and #2.

Although most of the structure of the two-particle vertex is indeed centered around zero frequency, its precise extent strongly depends on the cutoff scale $\Lambda$ (see, e.g., Figs. 4 and 8) and a static frequency grid will therefore fail to faithfully resolve the evolution of frequency structures under the fRG flow. Furthermore, multipeak structures that are present in several vertex components will in general not be captured by logarithmic sampling.

Fig. 9

Structure factor for the paramagnetic setup along a high-symmetry path of the cubic lattice Brillouin zone. The results are in good agreement between both codes, both for $\ell \!=\! 1$ and $\ell \!=\! 3$, showing that correlations are strongest around the ${\varvec{M}}$ point. Here, the peak sharpens with increasing loop order, and the $\ell =3$ flow predicts enhanced long-range correlations

Fig. 10

Decomposition of the $\smash {\gamma _{t, \langle ij\rangle }^{s}(\omega , \nu , \nu ')}$ vertex in the paramagnetic setup as in Fig. 5, for the $\ell = 3$ flows at $\Lambda / J = 0.3$. Here, all asymptotic classes are of the same order of magnitude, and structures with multiple peaks are present along the fermionic frequency cut (second column)

To address both shortcomings, codes #1 and #2 introduce hybrid frequency meshes using linear spacing around zero frequency augmented by an algebraic (code #1) or logarithmic (code #2) part to capture the high-frequency behavior in the asymptotic classes $K_n$ or $Q_n$. The parameters of these meshes are then independently rescaled for different vertex components making use of sophisticated scanning routines (see [1, 2] for further details).

4.2 Evaluation of bubble integrals

Having fixed the frequency discretization, the evaluation of frequency integrals in loop and bubble functions necessitates the use of a quadrature rule. In earlier implementations, a trapezoidal quadrature was used, with integration points coinciding with the frequency mesh of the vertex. As discussed above, this procedure yields good resolution around the origin of the integration variable. For $1\ell$ calculations, the bubble function consists of a single-scale and a full propagator, the former being more strongly peaked than the latter. As the integration variable was usually shifted such that the origin coincided with the more important pole of the single-scale propagator, at least the dominant contribution was accounted for in previous implementations.

Table 1 Technical summary of the algorithmic choices in code #1 and #2

In higher loops, however, both propagators enter the bubble on equal footing, necessitating adaptive routines to deal with the enriched frequency structure. This is illustrated in Fig. 11, where we compare the results of integrating the bare susceptibility

$$\begin{aligned} \chi ^{\Lambda }_{0}(\omega ) = \frac{1}{4 \pi } \int d\nu \, G^{\Lambda }_{0}(\nu + \tfrac{\omega }{2}) \, G^{\Lambda }_{0}(\nu - \tfrac{\omega }{2}) \,, \end{aligned}$$

Fig. 11

Evaluation of bubble integrals. Comparison of the bare susceptibility $\chi ^{\Lambda }_{0}(\omega ) = \frac{1}{4 \pi } \int d\nu \,G^{\Lambda }_{0}(\nu + \frac{\omega }{2})\, G^{\Lambda }_{0}(\nu - \frac{\omega }{2})$ obtained numerically via adaptive and static quadrature. The adaptive method utilizes the Simpson rule, while the static method applies a trapezoidal rule to a fixed logarithmic frequency discretization (see main text for more details). For frequencies larger than the scale set by the cutoff $\Lambda$, the non-adaptive integration becomes unstable and is plagued by rapid oscillations. By contrast, the adaptive routine yields stable results even beyond the small frequency regime and is therefore crucial to obtain accurate results for the vertex functions and their asymptotic behavior

i.e., the simplest bubble-like integral encountered during the fRG flow. Using trapezoidal quadrature over a fixed set of 60 logarithmically distributed integration points between $\nu _{\text {min}} = 10^{-3} J$ and $\nu _{\text {max}} = 250 J$, we find strong deviations for frequencies $\omega / \Lambda > rsim 1 \sim 10$ compared to the results produced with the adaptive routine of code #2 (see Ref. [2] for further details). Moreover, at small cutoffs $\Lambda / J \lesssim 1$, the non-adaptive result is plagued by rapid oscillations, rendering it numerically unstable and thus inapplicable. Analytically, an asymptotic falloff with a power law $\omega ^{-2}$ is expected, and this is reproduced perfectly by the adaptive integrator.

We emphasize that the test case considered here merely constitutes the simplest version of a bubble-like integral computed within the pffRG flow. In general, the propagators in bubble functions are dressed with self-energy insertions and additionally contracted with two-frequency dependent vertices. One should, therefore, expect even larger numerical errors for full fRG calculations that utilize non-adaptive quadrature.

Fig. 12

Scaling of relative runtime with numerical parameters. Median computational runtime of 60 samples of a single calculation of the right-hand side of the flow equation for $\Lambda / J = 1$ relative to the runtime of the fastest computation in each series. Calculations start from a parquet solution to make the code integrate over non-trivial frequency structures. The numerical parameters for all plots are fixed to $N_\omega = 50$, $N_\nu = 30$, $\xi = 4$ and $\ell = 1$, if not varied. The asymptotic behavior expected analytically is achieved in all cases (dashed red lines)

4.3 Flow integration

The integration of the RG flow can, in principle, be performed using any standard solver for ordinary differential equations. While earlier works used an Euler scheme with decreasing step-sizes (see, e.g., Ref. [46]), we employ higher order solvers in the Runge–Kutta family with adaptive step-size control to achieve maximum accuracy while being numerically efficient to operate. It is of particular importance to implement an error-controlling method near ordering instabilities such as the ferromagnetic setup in Sect. 3, as otherwise numerical errors may become unacceptably large even at scales $\Lambda \approx J$.

4.4 Initial condition

The final ingredient to set up the pffRG flow is an appropriate initial condition. In the UV limit $\Lambda \rightarrow \infty$, the pseudofermion vertex is given by the bare spin coupling, which, in numerical calculations, is naturally implemented using J as the initial condition at a large but finite value of $\Lambda$. The mfRG flow will, by construction, reproduce a solution to the parquet equations [33,34,35], given an initial condition consistent with them. Therefore, we solve the regularized parquet equations iteratively for an initial scale $\Lambda / J = 5$ and use the resulting self-energy and reducible vertices as a dynamic, i.e., frequency-dependent starting point for the fRG flow [1].

4.5 Scaling analysis

Most of the runtime needed to evaluate the right-hand side of the flow equations is spent calculating the derivative of the high-dimensional two-particle vertex as given in Eq. (6). In comparison, the computation time spent for the self-energy derivative of Eq. (5) is negligible. Consequently, the (asymptotic) computational complexity is given by

$$\begin{aligned} \mathcal {O}\left( N_\xi ^2 \times N_{\omega } N_\nu ^2 \times \ell \right) \,, \end{aligned}$$

where $N_\xi$ is the number of (symmetry reduced [1, 2]) lattice sites, $N_\omega$ ($N_\nu$) the number of bosonic (fermionic) frequencies, and $\ell$ denotes the number of loops. The total number of sites, in turn, is expected to follow a $\mathcal {O}(\xi ^d)$ dependence, where $\xi$ is the maximal correlation length considered and d is the spatial dimensionality of the underlying lattice, with $d=3$ for the simple cubic lattice at hand.

To demonstrate that we indeed reach this asymptotic algorithmic scaling also in numerical implementations we show, in Fig. 12, the median runtime data for 60 evaluations of the right-hand side of the fRG equations obtained using code #2. For the number of bosonic and fermionic frequencies, the expected linear and quadratic behavior, respectively, is achieved over the whole parameter range. Note that, due to the adaptive integration and parallelization used, slight deviations from the theoretical scaling are to be expected. Similarly, the scaling in the maximal correlation length $\xi$ is achieved for the whole parameter range. In the number of loops, the linear scaling sets in at $\ell =5$, while for smaller $\ell$ a steeper slope is found. We attribute this behavior to the contributions of higher loops becoming successively smaller, leading to faster converging adaptive loop integrals for given absolute and relative tolerances. That way, the initial overhead of computing two (three) loop corrections, which require twice (thrice) the number of integrals to be evaluated compared to $\ell = 1$, diminishes with increasing loop number and the analytically expected scaling, linear in $\ell$, is recovered.

Table 2 Number of (symmetry reduced) vertex flow equations for Heisenberg models on the cubic lattice as a function of the maximum correlation length $\xi$. The number of positive frequencies is fixed to 60 (50) for the bosonic (fermionic) Matsubara axis

As a final remark, we mention that the number of vertex flow equations, another measure of algorithmic complexity, grows rapidly as one increases the maximal correlation length considered for a given lattice model. This is summarized in Table 2.

5 Conclusions

We benchmarked two state-of-the-art codes for solving pseudofermion functional renormalization group equations. Our analysis considered both physical observables, i.e. spin-spin correlation functions and structure factors, as well as fermionic vertex functions (self-energy and two-particle vertex) for ferro- and antiferromagnetic models on the simple cubic lattice.

For the nearest-neighbor ferromagnet, both codes were in quantitative agreement at least until $\Lambda / J > rsim 0.76$, where they consistently predicted a breakdown of the RG flow, indicated by a sharp peak (for $\ell = 1$) or a divergence (for $\ell = 3$) in the spin-spin correlations. The energy scale $\Lambda _{\mathrm {T}}$ associated with this numerical instability slightly differed, which necessitated an in-depth comparison of the influence of the numerical frequency grid on the obtained results. We found that both fRG solvers, due to the emergence of a singular peak in the t reducible vertex functions, become sensitive to the precise mesh spacing and thus predict marginally different critical scales, although the physical conclusion drawn from the RG flow, i.e. the onset of long-range ferromagnetic order, remains the same.

For the antiferromagnetic setup, the $\ell = 1$ results obtained by both codes were in agreement with one another and previous studies [21], predicting a paramagnetic state, signified by a regular RG flow down to the infrared. For $\ell = 3$, similar numerical agreement between the two codes was found. However, the physical results changed qualitatively: the flow of the spin-spin correlator diverged around $\Lambda / J \approx 0.24$, accompanied by sharp Bragg peaks at the ${\varvec{M}}$ points indicating the formation of antiferromagnetic order at low temperatures. This reinstantiates the importance of including higher loop corrections in pffRG to avoid overestimating the extent of paramagnetic phases and to obtain more accurate predictions of ground states in frustrated quantum magnets.

We also elaborated on the importance of employing adaptive numerical algorithms to obtain robust results at all stages of the flow. More explicitly, there are extended structures with multiple peaks in the three-dimensional frequency dependence of several vertex components. As these structures are sizable, it is crucial to resolve them in an accurate manner. We found fixed logarithmic frequencies to be insufficient for structures not centered at zero frequency, and rely instead on adaptive frequency meshes that have been specifically optimized for pffRG vertices. Furthermore, we demonstrated that the commonly employed quadrature of a trapezoidal rule over a static, logarithmic mesh fails to produce the analytically expected behavior of bare bubble integrations at large frequencies. It is thus unsuitable for providing the essential Matsubara integrals for error-controlled fRG flows. By contrast, the implementations presented and benchmarked here solve these problems using highly accurate, yet efficient adaptive routines (see Table 1). We thus believe that, moving forward, they will be widely used for unbiased calculations of (multiloop) ground-state phase diagrams of frustrated magnets from pffRG.