Cluster functional renormalization group and absence of a bilinear spin liquid in the $J_1$-$J_2$-Heisenberg model

TL;DR

This paper introduces a cluster pf-FRG method to better identify spin liquid states in frustrated quantum magnets, applying it to the $J_1$-$J_2$ Heisenberg model and finding no evidence for a bilinear spin liquid in this system.

Contribution

The paper develops a cluster pf-FRG approach that links microscopic models to specific spin liquid states, improving the detection of such states compared to previous methods.

Findings

Correctly captures the $ ext{pi}$-flux spin liquid in the large-$N$ limit.

Suggests the $J_1$-$J_2$ Heisenberg model does not support fermion bilinear spin liquids.

Provides a detailed phase diagram near the frustrated regime for $N=2$.

Abstract

The pseudofermion functional renormalization group (pf-FRG) has been put forward as a semi-analytical scheme that, for a given microscopic spin model, allows to discriminate whether the low-temperature states exhibit magnetic ordering or a tendency towards the formation of quantum spin liquids. However, the precise nature of the putative spin liquid ground state has remained hard to infer from the original (single-site) pf-FRG scheme. Here we introduce a cluster pf-FRG approach, which allows for a more stringent connection between a microscopic spin model and its low-temperature spin liquid ground states. In particular, it allows to calculate spatially structured fermion bilinear expectation values on spatial clusters, which are formed by splitting the original lattice into several sublattices, thereby allowing for the positive identification of a family of bilinear spin liquid states. As an application of this cluster pf-FRG approach, we consider the $J_1$-$J_2$ SU($N$)-Heisenberg model on a square lattice, which is a paradigmatic example for a frustrated quantum magnet exhibiting quantum spin liquid behavior for intermediate coupling strengths. In the well-established large-$N$ limit of this model, we show that our approach correctly captures the emergence of the $\pi$-flux spin liquid state at low temperatures. For small $N$, where the precise nature of the ground state remains controversial, we focus on the widely studied case of $N=2$, for which we determine the low-temperature phase diagram near the strongly-frustrated regime after implementing the fermion number constraint by the flowing Popov-Fedotov method. Our results suggest that the $J_1$-$J_2$-Heisenberg model does not support the formation of a fermion bilinear spin liquid state.

Figures

Tables

Table 1. Table 1: Momentum structure and non-vanishing entries of the sublattice-space matrices η X / Y superscript 𝜂 𝑋 𝑌 \eta^{X/Y} for the initial interaction of the J 1 subscript 𝐽 1 J_{1} - J 2 subscript 𝐽 2 J_{2} model.

Coupling	${\eta }^X$	${\eta }^Y$	$f(p_1,p_2,p_3,p_4)$
$J_{1,\Lambda }$	13	31	$\cos \left(p_{2,x}-p_{3,x}\right)$
$J_{1,\Lambda }$	14	41	$\cos \left(p_{2,y}-p_{3,y}\right)$
$J_{1,\Lambda }$	23	32	$\cos \left(p_{2,y}-p_{3,y}\right)$
$J_{1,\Lambda }$	24	42	$\cos \left(p_{2,x}-p_{3,x}\right)$
$J_{2,\Lambda }$	12	21	$\cos \left(p_{2,x}-p_{3,x}\right)\cos \left(p_{2,y}-p_{3,y}\right)$
$J_{2,\Lambda }$	34	43	$\cos \left(p_{2,x}-p_{3,x}\right)\cos \left(p_{2,y}-p_{3,y}\right)$

Table 2. Table 2: Momentum structure and non-vanishing entries of the sublattice-space matrix η X superscript 𝜂 𝑋 \eta^{X} for density-type order parameters.

Order Parameter	${\eta }^X$	$f\left(p\right)$
$Q^{\left(1\right)}$	{13,31}	$\cos \left(p_x\right)$
$Q^{\left(2\right)}$	{14,41}	$\cos \left(p_y\right)$
$Q^{\left(3\right)}$	{23,32}	$\cos \left(p_y\right)$
$Q^{\left(4\right)}$	{24,42}	$\cos \left(p_x\right)$
$Q^{\left(5\right)}$	{12,21}	$\cos \left(p_x\right)\cos \left(p_y\right)$
$Q^{\left(6\right)}$	{34,43}	$\cos \left(p_x\right)\cos \left(p_y\right)$

Table 3. Table 3: Full interaction structure of the symmetric two-sublattice J 1 subscript 𝐽 1 J_{1} model. Matrix definitions are: σ ± = 1 2 ​ ( σ x ± σ y ) superscript 𝜎 plus-or-minus 1 2 plus-or-minus superscript 𝜎 𝑥 superscript 𝜎 𝑦 \sigma^{\pm}=\frac{1}{2}(\sigma^{x}\pm\sigma^{y}) and σ o ​ u = 1 2 ​ ( 𝟙 ± σ z ) superscript 𝜎 𝑜 𝑢 1 2 plus-or-minus 1 superscript 𝜎 𝑧 \sigma^{ou}=\frac{1}{2}(\mathbbm{1}\pm\sigma^{z}) with the usual definitions of the Pauli matrices σ x / y / z superscript 𝜎 𝑥 𝑦 𝑧 \sigma^{x/y/z} .

Coupling	${\eta }^X$	${\eta }^Y$	$f(p_1,p_2,p_3,p_4)$
$J_{\Lambda }$	${\sigma }^{+}$	${\sigma }^{-}$	$\cos \left(p_{2,x}-p_{3,x}\right)+\cos \left(p_{2,y}-p_{3,y}\right)$
$L_{\Lambda }$	${\sigma }^o$	${\sigma }^u$	$\cos \left(p_{3,x}-p_{4,x}\right)+\cos \left(p_{3,y}-p_{4,y}\right)$
$X_{\Lambda }$	${\sigma }^o$	${\sigma }^o$	1
$X_{\Lambda }$	${\sigma }^u$	${\sigma }^u$	1

Table 4. Table 4: New interaction structures due to finite pairing order parameters. Spin space is represented by the diagonal matrix δ 𝛿 \delta and the totally antisymmetric matrix ϵ italic-ϵ \epsilon . The geometric structure f k subscript 𝑓 𝑘 f_{k} is given in Eq. ( 20 ).

Type	Grassmann Structure	Momentum Structure
Original	$\left({\Psi }^{†}{\eta }^X\delta \Psi \right)\left({\Psi }^{†}{\eta }^Y\delta \Psi \right)$	$f_k(1,-1,1)$
Mixed	$\left({\Psi }^{†}{\eta }^Xϵ{\Psi }^{\ast }\right)\left({\Psi }^{†}{\eta }^Y\delta \Psi \right)$	$f_k(-1,-1,1)$
	$\left({\Psi }^T{\eta }^Xϵ\Psi \right)\left({\Psi }^{†}{\eta }^Y\delta \Psi \right)$	$f_k(1,-1,1)$
Pairing	$\left({\Psi }^{†}{\eta }^Xϵ{\Psi }^{\ast }\right)\left({\Psi }^{†}{\eta }^Yϵ{\Psi }^{\ast }\right)$	$f_k(1,1,1)$
	$\left({\Psi }^T{\eta }^Xϵ\Psi \right)\left({\Psi }^T{\eta }^Yϵ\Psi \right)$	$f_k(1,1,1)$
	$\left({\Psi }^{†}{\eta }^Xϵ{\Psi }^{\ast }\right)\left({\Psi }^T{\eta }^Yϵ\Psi \right)$	$f_k(1,-1,-1)$

Equations

$$\mathcal{H}=\frac{J_1}{N}\sum _{⟨ij⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+\frac{J_2}{N}\sum _{⟨⟨ij⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$${\partial }_{\Lambda }{\Gamma }_{\Lambda }=\frac{1}{2}{\int }_{\tau }\sum _i{f}^{†}{\partial }_{\Lambda }{\mathcal{P}}_{\Lambda }f+\frac{1}{2}\mathrm{STr}\left[\frac{{\partial }_{\Lambda }{\mathcal{P}}_{\Lambda }}{{\Gamma }_{\Lambda }^{(2)}}\right],$$

$$S_i^{\mu }=f_{i\alpha }^{†}{T}_{\alpha \beta }^{\mu }{f}_{i\beta }^{},$$

$$f_{i\alpha }^{†}{f}_{i\alpha }^{}=\frac{N}{2}\forall i.$$

$$S$$

$$-\frac{J_1}{N}\sum _{⟨ij⟩}{f}_{i\alpha }^{†}{f}_{j\alpha }^{}{f}_{j\beta }^{†}{f}_{i\beta }^{}-\frac{J_2}{N}\sum _{⟨⟨ij⟩⟩}{f}_{i\alpha }^{†}{f}_{j\alpha }^{}{f}_{j\beta }^{†}{f}_{i\beta }^{}],$$

$${\Psi }^{†}=\left(f^{†,A},f^{†,B}\right)\to {\Psi }^{†}=\left(f^{†,A},f^{†,B},f^{†,C},f^{†,D}\right).$$

$$\frac{J_n}{N}{\int }_{p_1..p_4\in \mathrm{TZ}}\left({\Psi }_{p_1}^{†}{\eta }_n^X{\Psi }_{p_2}^{}\right)\left({\Psi }_{p_3}^{†}{\eta }_n^Y{\Psi }_{p_4}^{}\right)$$

$$\cdot f(p_1,p_2,p_3,p_4)\cdot {\delta }_{p_1{p}_2{p}_3{p}_4},$$

$$Q^{(n)}{\int }_{p\in \mathrm{TZ}}\left({\Psi }_p^{†}{\eta }_n^X{\Psi }_p^{}\right)f(p).$$

$${\mathrm{Gen}}_{\mathrm{SU}(4)}=\{{\hat{\lambda }}_1,{\hat{\lambda }}_6,{\hat{\lambda }}_7,{\hat{\lambda }}_{12},{\hat{\lambda }}_{13},{\hat{\lambda }}_{14}\}.$$

$${\hat{Q}}_{\Lambda \to {\Lambda }_{\mathrm{UV}}}^{\mathrm{BZA}}=Q_{\mathrm{UV}}\cdot \left(\begin{array}{cc}0& \text{mathbbm}{1}_2+{\sigma }_x\\ \text{mathbbm}{1}_2+{\sigma }_x& 0\end{array}\right).$$

$${\hat{\lambda }}_{\pi }={\hat{\lambda }}_1+{\hat{\lambda }}_6,$$

$${\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}R}$$

$${\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}I}$$

$${\hat{Q}}_{{\Lambda }_{\mathrm{UV}}}^{\mathrm{BZA}}+\sum _i{Q}_i{\hat{\lambda }}_i,{\hat{\lambda }}_i\in \left\{{\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}R}\cup {\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}I}\right\},$$

$${\Psi }^{†}{\hat{\lambda }}_i\Psi ,{\hat{\lambda }}_i\in \left\{{\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}R}\cup {\mathrm{Gen}}_{\mathrm{flux},\text{mathbbm}I}\right\},$$

$${\mu }_{\mathrm{PF}}=i\frac{\pi T}{2},$$

$${\Gamma }_{\Lambda }^{(2\Psi )}=\sum _n{\int }_{p\in \mathrm{TZ}}\left[{\Psi }_p^{†}\left(i{\omega }_n{\vartheta }_{\Lambda }^{-1}+{\mu }_{\Lambda }\right)\text{mathbbm}1{\Psi }_p\right],$$

$${\vartheta }_{\Lambda }=\{\begin{array}{cc}0& |{\omega }_n|<\Lambda -\pi T\\ \frac{1}{2}+\frac{|{\omega }_n|-\Lambda }{2\pi T}& \Lambda -\pi T\le |{\omega }_n|\le \Lambda +\pi T\\ 1& \Lambda +\pi T\le |{\omega }_n|\end{array}.$$

$$S_i^{\mu }{S}_j^{\mu }=-\frac{1}{2}{f}_{i\alpha }^{†}{f}_{j\alpha }^{}{f}_{j\beta }^{†}{f}_{i\beta }^{}-\frac{1}{4}{f}_{i\alpha }^{†}{f}_{i\alpha }^{}{f}_{j\beta }^{†}{f}_{j\beta }^{},$$

$${\Gamma }_{\Lambda }^{(2\Psi )}=$$

$$+\sum _{i=1}^6{\Psi }_{p\alpha }^{†}\left(Q_{\Lambda }^{(i)}{\eta }_i+Q_{\Lambda }^{(i),\ast }{\eta }_i^T\right){\delta }_{\alpha \beta }\mathfrak{Q}(p){\Psi }_{p\beta }$$

$$+\sum _{j=1}^6{\Delta }^{(j),\ast }{\Psi }_{p\alpha }^{†}{\eta }_j^T{\Psi }_{-p\beta }^{\ast }{ϵ}_{\alpha \beta }\mathfrak{D}(p)$$

$$\displaystyle+\sum_{j=1}^{6}\Delta^{(j)}\Psi_{p\alpha}^{T}\eta_{j}\Psi_{-p\beta}\epsilon_{\alpha\beta}\mathfrak{D}(p)\Big{]}\,.$$

$$f_k(i,j,l)=$$

$$\cdot \cos (i{a}_y{k}_2+j{b}_y{k}_3+l{c}_y{k}_4).$$

$${\Gamma }_{\Lambda }^{(\gamma )}=\sum _n{\int }_{p\in TZ}i\gamma {\Psi }_p^{†}\mathrm{sign}({\omega }_n)\text{mathbbm}1{\Psi }_p^{}$$

$${\mathcal{F}}_{\Lambda }=\frac{J_n}{N}{\delta }_{m,n}\delta (\mathbf{p}-\mathbf{q})\left(\begin{array}{cc}\mathfrak{A}& \mathfrak{B}\\ \mathfrak{C}& \mathfrak{D}\end{array}\right),$$

$$\mathfrak{A}=$$

$$-({\Psi }_{\gamma }^{†}{\eta }^Y{)}^T({\Psi }_{\delta }^{†}{\eta }^X)f(0,p,0,-p),$$

$$\mathfrak{B}=$$

$$+({\Psi }_{\gamma }^{†}{\eta }^Y{)}^T({\eta }^X{\Psi }_{\delta }{)}^{T}f(-p,0,0,-p)$$

$$-{{\eta }^X}^T{\delta }_{\gamma \delta }({\Psi }_{\alpha }^{†}{\eta }^Y{\Psi }_{\alpha })f(-p,-p,0,0)$$

$$-{{\eta }^B}^T{\delta }_{\gamma \delta }({\Psi }_{\alpha }^{†}{\eta }^X{\Psi }_{\alpha })f(0,0,-p,-p),$$

$$\mathfrak{C}=$$

$$+({\eta }^Y{\Psi }_{\gamma })({\Psi }_{\delta }^{†}{\eta }^X)f(0,p,p,0)$$

Full text

I Introduction

Quantum magnets are host to an astounding range of fascinating phenomena which go beyond the realm of classical magnetism. This includes the formation of multipolar order, such as spin nematic states, or more generally topological order, along with the appearance of fractionalized excitations and long-range entanglement Lacroix et al. (2011). The latter are hallmarks of quantum spin liquids (QSLs) Balents (2010), which have been widely studied in the context of frustrated magnets. Substantial progress in the understanding of the emergence of such QSLs has been achieved through the seminal work by Kitaev Kitaev (2006), providing an exactly solvable microscopic model with various QSL ground states. Aside from this particular model system, it remains, however, notoriously difficult to establish reliable connections between microscopic spin Hamiltonians and their possible QSL ground states.

A paradigmatic example for the challenges to be met is the $J_{1}$-$J_{2}$ Heisenberg model. Here, spins on a two-dimensional square lattice interact antiferromagnetically with their nearest and next-to-nearest neighbors with couplings $J_{1}$ and $J_{2}$, respectively, and the level of frustration can be tuned by the ratio of the couplings $J_{2}/J_{1}$. The apparent simplicity of this model combined with an early finding Chandra and Doucot (1988) of a low-temperature non-magnetic phase in the strongly frustrated intermediate coupling regime has drawn the attention of researchers for over three decades now. Many methods have been developed and applied to clarify the precise nature of the frustrated ground state, including exact diagonalization Einarsson and Schulz (1995); Schulz et al. (1996); Capriotti and Sorella (2000); Mambrini et al. (2006); Richter and Schulenburg (2009), density matrix renormalization group Jiang et al. (2012); Gong et al. (2014); Wang and Sandvik (2018), tensor network techniques Wang et al. (2016); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Liu et al. (2018), variational approaches Mezzacapo (2012); Hu et al. (2013), different expansion schemes Chandra and Doucot (1988); Sirker et al. (2006); Aktersky and Syromyatnikov (2016) and others Beach (2009); Ren et al. (2014); Metavitsiadis et al. (2014); Chou and Chen (2014); Richter et al. (2015). For the ground state, gapped Jiang et al. (2012) and gapless Hu et al. (2013); Richter et al. (2015); Chou and Chen (2014); Wang et al. (2016); Wang and Sandvik (2018) spin liquids have been found as well as different types of valence bond solids Capriotti and Sorella (2000); Mambrini et al. (2006); Sirker et al. (2006); Beach (2009); Gong et al. (2014); Metavitsiadis et al. (2014); Aktersky and Syromyatnikov (2016); Haghshenas and Sheng (2018), but the results are widespread and not congruent within or across the applied approaches, leaving a rather unsatisfactory situation.

The large range of incompatible results on the ground state of the frustrated $J_{1}$-$J_{2}$ Heisenberg model suggests that it is rather sensitive to any bias introduced by an approximate approach. The goal of this work is to avoid any such bias by construction. We employ and expand renormalization group techniques, which have successfully been used in diverse physical contexts to systematically study effective low-energy models of microscopic theories. A particular realization, the pseudofermion functional renormalization group (pf-FRG) has the capacity to deal with strongly correlated spin systems and many works have proven its ability to accurately characterize quantum magnets in complicated lattice geometries with respect to ordering patterns as well as critical temperatures and couplings Reuther and Wölfle (2010); Reuther and Thomale (2011); Reuther et al. (2011); Reuther and Thomale (2014); Buessen and Trebst (2016); Iqbal et al. (2016); Baez and Reuther (2017); Buessen et al. (2018a).

Until recently, a central drawback of the pf-FRG approach has been that the occurrence of magnetized and even certain types of QSL ground states leads to an instability in the renormalization group flow. While this instability itself represents a hallmark of ordering, it prevents a continuation of the flow towards the low-energy regime. Therefore, the precise nature of the ground state or even competing orders remained elusive. In Refs. Buessen et al., 2018b; Roscher et al., 2018 it was demonstrated, that the pf-FRG can be properly regularized such that the occurrence of an instability does not lead to a breakdown of the description itself anymore. Since this result can be achieved with arbitrarily small (initial) bias, it can be used to positively identify and characterize spin-liquid phases that can be represented by expectation values of pseudofermion bilinears Wen (2002), which we refer to as bilinearly ordered spin liquids in the following.

In this work, we substantially extend this regularized pf-FRG approach by working with multi-site clusters in lieu of the single-site perspective employed hitherto. Importantly, this allows to map part of the spatial structure of spin-liquid order parameters onto the pseudofermion representation itself. This, in turn, enables us to systematically include more spin-liquid ordering patterns in the characterization of the low-energy physics. Employing this cluster pf-FRG approach to the $J_{1}$-$J_{2}$ SU($N$)-Heisenberg model on the square lattice, we show that benchmark calculations for the large-$N$ case indeed exhibit the expected emergence of a $\pi$-flux spin liquid state in the low-energy regime. This reveals the true power of the cluster pf-FRG approach: the RG flow automatically and reliably chooses the energetically most favorable ordering pattern available at the given level of approximation. This distinct feature of the approach can not only be used to systematically improve the quality of predictions for spin-liquid phases, but statements about their very existence are found to be independent of the level of approximation employed for the spatial structure. Taking into account the full set of non-magnetic bilinear ordering patterns, the cluster pf-FRG approach thus allows us to obtain definite statements on the existence of magnetically ordered or spin-liquid phases in the finite-temperature phase diagram of the $J_{1}$-$J_{2}$ model.

I.1 Line of Arguments and Overview of Results

We develop a generalization of the pf-FRG method to the $J_{1}$-$J_{2}$ Heisenberg model at zero and finite temperatures, aiming at a clarification of whether the frustrated regime can be described by a bilinear spin liquid. First, in order to establish a well-controlled limit without having to deal with complicated fermion number constraints, we revisit the SU($N$) symmetric model at large $N$. We show that our generalization correctly captures the emergence of the $\pi$-flux spin liquid state at low energies, cf. Sec. III, by introducing a splitting of the square lattice into four sublattices, cf. Sec. II.2, and an appropriate set of order parameters. In particular, we show that the RG flow automatically selects the energetically favored ground state. This analysis corroborates the suitability of our approach to detect any (non-magnetic) bilinearly ordered spin liquid states despite the employed approximations. Moreover, we clearly exhibit that the spin-liquid behavior at large $N$ results from the renormalization group flow of the four-fermion function and is not related to a frequency-dependent self-energy, Sec. IV.4. Introducing an artificial damping into the self-energy causes an unphysical shift of the critical temperature and order parameter. With this knowledge, we turn to the physically relevant case of $N=2$, where we first carefully implement the fermion number constraint in terms of a flowing Popov-Fedotov method, Sec. IV.1, and then calculate the low-temperature phase diagram as a function of the ratio $J_{2}/J_{1}$, Sec. IV. We find that our approach predicts the appearance of magnetic instabilities through the full range of $J_{2}/J_{1}$ with only a moderate suppression in the strongly-frustrated regime. In view of the insight on the role of frequency dependencies, this suggests that the $J_{1}$-$J_{2}$ Heisenberg model does not support the formation of a bilinearly ordered spin liquid state. Conclusions and future prospects are given in Sec. V.

II Model and Implementation

The goal of our analysis is to characterize the ground state of the Heisenberg-Hamiltonian

[TABLE]

implemented on a two-dimensional square lattice. The nearest ($J_{1}$) and next-to-nearest ($J_{2}$) neighbor couplings are positive and thus favor antiferromagnetic correlations. Their ratio $g=J_{2}/J_{1}$ is an important quantity, encoding the frustration of the system.

II.1 Pseudofermion FRG

Our tool of choice is the pseudofermion formulation of functional renormalization group (FRG). At its heart lies the functional differential Wetterich equation Wetterich (1993)

[TABLE]

for the running effective action $\Gamma_{\Lambda}$ with multiplicative regularization Salmhofer and Honerkamp (2001). Here, $\{f^{\dagger},f\}$ are anticommuting Grassmann numbers. At the initial scale $\Lambda\rightarrow\Lambda_{\mathrm{UV}}$, $\Gamma_{\Lambda}$ is given by the microscopic action corresponding to the Hamiltonian (1). After choosing a regulator function by which to multiply the free propagator $\mathcal{P}_{\Lambda}$, solving the flow equation, Eq. (2), in principle yields the full effective action $\Gamma_{\Lambda\rightarrow 0}$ from which thermodynamic properties can be read off.

Practically, two prerequisites are needed to actually perform this procedure. Firstly, the Hamiltonian operator (1) needs to be recast in a form that is amenable to our chosen formulation of the FRG flow Eq. (2). To that end, we choose the pseudofermion representation of spin operators Abrikosov (1965)

[TABLE]

where $T^{\mu}$ are matrices of the underlying spin symmetry group. Thus, $T^{\mu}=\frac{1}{2}\sigma^{\mu}$ for the physical SU(2) symmetry with $\sigma^{\mu}$ being Pauli matrices, and generalized Gell-Mann matrices for general SU($N$). To eliminate unphysical degrees of freedom, the (pseudo)fermions are also subject to the local constraint on the particle number,

[TABLE]

Making use of commutation relations between the $T^{\mu}$ and applying the constraint to eliminate particle number operators, the Hamiltonian (1) can be recast into an action that serves as initial condition for the flow equation,

[TABLE]

where, in the following, the explicit spin indices will often be suppressed for better readability.

We note that the flow equation (2) can usually not be solved without an approximate ansatz for $\Gamma_{\Lambda}$. The choice of this ansatz is crucial since all structures relevant to the physics of interest need to be included, while keeping the computational cost manageable. Systematic development and benchmarking of this truncation for $\Gamma_{\Lambda}$ is a central objective of our analysis in the forthcoming sections.

II.2 Sublattice representation: Spatially structured orders

The inclusion of explicit flowing order parameters in the pf-FRG scheme by construction reduces the symmetry of the problem and thus increases computational complexity. For this reason, a point-like momentum-space projection was applied in Ref. Roscher et al., 2018 instead of the conventional spatially resolved pf-FRG schemes Reuther and Wölfle (2010). In this projection scheme, $J_{1,2}$ does not depend on frequency or momentum. This approach trades spatio-temporal resolution for explicit information on ordering tendencies and access to the physically relevant low-energy regime. Consequently, only a spatially homogeneous Baskaran-Zou-Anderson (BZA)-phase could be detected in the large-$N$ analysis so far, while it is known Affleck and Marston (1988); Arovas and Auerbach (1988) that (spatially structured) $\pi$-flux or even Peierls phases are energetically preferred.

For the analysis of possible spin liquid phases in the $J_{1}$-$J_{2}$ model, this seems too restrictive, as it would exclude the majority of possible spin-liquid ordering patterns Wen (2002). Furthermore, the next-to-nearest neighbor interaction $\sim J_{2}$ needs to be included and discriminated from the $\sim J_{1}$ structure. In order to discriminate the $J_{1}$- and $J_{2}$-coupling structures or different spatially-structured ordering, one could abandon the simple point-like projection in momentum space in favor of an actual, finite discretization of the Brillouin zone. This would increase computational cost to a degree approaching or even exceeding that of the conventional pf-FRG scheme.

We develop a different approach here, where the point-like projection scheme is kept and discrimination between structures on different bonds can be performed algebraically. When representing the nearest-neighbor structure of the $J_{1}$ interaction, it is natural to artificially split the original square lattice into two equivalent sublattices $A$ and $B$, defining the interaction on the links between them Arovas and Auerbach (1988). We perpetuate this idea by moving to four artificial sublattices, see Fig. 1 . Spinors can now be defined on the corresponding four-dimensional space:

[TABLE]

By properly choosing $4\times 4$ matrices $\eta^{X/Y}$, the general interaction structure in frequency-momentum space,

[TABLE]

can now be restricted to nearest- or next-to-nearest neighbors. Note that we have suppressed Matsubara labels. Further, TZ defines the “tiny zone”, corresponding to one quarter of the Brillouin zone of the original lattice. The geometry-induced momentum structure is represented by $f(p_{1},p_{2},p_{3},p_{4})$ and $\delta_{p_{1}p_{2}p_{3}p_{4}}$ is a momentum-conserving $\delta$-function.

Using the representation in Eq. (7) it is possible to resolve spatial patterns of intermediate range in a simple algebraic and systematic way without having to consider the explicit momentum dependence of interaction vertices and self-energies. One can, in principle, always increase the number of artificial sublattices to increase spatial range. Furthermore, this strategy can be generalized to other lattice geometries. For example, a 2D Kagome lattice may be described as a triangular lattice, split into four sublattices, one of which is eliminated Betts (1995). This will be the subject of future work. We also note that the present scheme allows for a comparatively simple algebraic analysis of spatial symmetry breaking patterns. This aspect will be discussed in detail in Sec. III.2, below.

The interaction structure given by the Hamiltonian (1) with respect to the sublattice definition from Eq. (6) can now be described uniquely by the non-vanishing entries of the respective $\eta^{X/Y}$ matrices, see Tab. 1.

The geometric momentum structure follows directly from the spatial information contained in the $\eta^{X/Y}$ matrices. Here, we exclusively consider couplings between the nearest available sublattice sites for a given $\eta^{X/Y}$-matrix. This amounts to an approximation, which will be discussed in Sec. IV.2, since for finite $N$, longer-ranged contributions to the respective vertices are generated. We disregard those to keep the system of RG equations closed.

III Large-$N$ analysis

As a first step towards a systematic study of the full model, Eq. (1), we perform a large-$N$ analysis. This provides useful hints on the character of a possible spin liquid phase in the frustrated regime and guides the more intricate search at the physical value of $N=2$. Furthermore, for large $N$, it is sufficient to implement the fermion number constraint (4) on average only Arovas and Auerbach (1988), which greatly simplifies the analysis.

We consider the well-understood limit $J_{2}=0$ first: To regularize divergences in the RG flow of the running coupling $J_{1,\Lambda}$, we introduce bilinear density-type symmetry-breaking order parameters $Q^{(1..4)}$. These may now be addressed in a manner analogous to Eq. (7),

[TABLE]

Possible order parameters compatible with four sublattices are given by their matrix structure in Tab. 2.

At large-$N$, pairing-type order parameters $\Delta$ are suppressed and will therefore not be considered here. When we proceed to $N=2$ below, we will take them into account, though. Moreover, since we are primarily interested in identifying spin-liquid states, we will not introduce magnetic order parameters $\sim M^{(n)}\Psi^{\dagger}_{\alpha}\tau_{\alpha\beta}\eta^{X}\Psi_{\beta}$ where $\tau_{\alpha\beta}\neq\delta_{\alpha\beta}$ in spin space, neither for $N=2$ nor for large $N$. This encompasses all possible bilinear order parameters. Magnetic phases will therefore be signaled by a divergence of the RG flow. In contrast, bilinearly-ordered spin-liquid states are identified by a regular RG flow with a finite order parameter in the infrared.

Due to the presence of the symmetry-breaking order parameters, new interaction structures besides the initial ones in Tab. 1 are generated during the flow. In the space of four sublattices, up to $4^{2}\cdot 4^{2}=256$ different structures can in principle occur and may be discriminated algebraically. However, not all of these are finite in the present setup. It can be shown that only matrices $\eta^{X/Y}$ which are already present in the initial $\Gamma_{\Lambda\rightarrow\Lambda_{\mathrm{UV}}}$ are available for combination among themselves at $N\rightarrow\infty$, cf. App. A for details. This reduces the number of newly generated interaction structures down to 36. Some of theses are related by hermitean conjugation, which can be used to further reduce the number of differential equations. The resulting 40 coupled differential equations are rather complex. We therefore refrain from writing them down explicitly. Determining potentially non-vanishing couplings as well as the construction of the flow equations themselves follows a set of comparatively simple and straightforward rules. Both problems are therefore amenable to automated and/or numerical analysis.

III.1 Staggered-flux spin liquid ground state

For large $N$, the mean-field approach to the SU($N$) Heisenberg model becomes exact and a staggered- or $\pi$-flux spin liquid phase arises as the ground state Affleck and Marston (1988); Wen (2002). The large-$N$ FRG flow equations are known to exactly reproduce the mean-field results if the Katanin scheme Katanin (2004) is employed, as shown previously Salmhofer et al. (2004); Roscher et al. (2018). Initializing the $Q^{(1..4)}$ with small $(\mathcal{O}(10^{-4}))$, generally complex, random numbers, we find excellent agreement for the absolute value of the order parameters as well as the magnetic susceptibility $\chi_{\mathrm{mag}}$ with the mean-field solution, see Fig. 2. The specific values of the individual order parameters themselves do depend on the initial values, see Fig. 3 for an example, in the sense that different but equivalent configurations arise: Except for one special configuration, see Sec. III.2, we always find an overall phase of $\pi$ around each plaquette as expected. The RG flow thus automatically selects the appropriate state.

We conclude that the sublattice splitting is indeed a suitable method to systematically capture spatial structures of spin liquid phases. A splitting into eight sublattices should give access to the true large-$N$ ground state which is known to be a Peierls phase Affleck and Marston (1988). We refrain from pursuing this, here, and focus on the $N=2$ case, in the following.

III.2 Algebraic ordering analysis

Xiao-Gang Wen’s symmetry-based classification of spin liquid states relies on a decoupling of fermion bilinears in terms of order-parameter-like mean-fields Wen (2002) which could, in principle, be related to the fermionic interaction by a Hubbard-Stratonovich transformation. In our fermionic RG formulation of spin models, such ordering tendencies manifest themselves in a divergence of the fermionic couplings Braun (2012); Buessen et al. (2018b) – unless regularized. This divergence is crucial for our approach, as it guarantees that this type of bilinear ordering cannot be missed. However, an order parameter which regularizes the fermionic RG flow, allowing for its continuation towards the infrared, may still not reflect the true symmetry-breaking pattern of the ground state. For instance, including the homogeneous BZA state in a simple two-sublattice setup leads to a convergent RG flow and no sign of something being amiss could be discerned. On the other hand, the ground state of the system is better described by the $\pi$-flux phase, as described in the previous section, which is characterized by an additional breaking of translation invariance. The difference between the two states lies in the spatial structure of the bilinear expectation value $Q_{ij}\sim\langle f^{\dagger}_{i}f^{\phantom{\dagger}}_{j}\rangle$, serving as an order parameter for translational/rotational symmetry breaking. The occurrence of this spatial structure by itself is not a priori related to a diverging coupling and can thus be missed by a naïve implementation. Considering the importance of additional breaking of discrete as well as non-compact symmetries for the classification of spin liquid states Wen (2002), this is an issue that will therefore be addressed in the following.

Employing the cluster pf-FRG introduced above, this situation is remedied in a systematic way. The square-lattice Heisenberg model is symmetric under primitive translations with respect to the lattice spacing $a$ along $\vec{e}_{x}$ and $\vec{e}_{y}$ and a number of other symmetries, forming a discrete subgroup of the non-compact Galilei group associated with non-relativistic continuum systems. When introducing a number $M$ of sublattices, part of this group is mapped onto a subgroup of the compact SU($M$)111As long as no pairing order parameters $\Delta\sim f_{i\alpha}f_{i\beta}\epsilon_{\alpha\beta}$ are introduced, the action may even be invariant under U($M$). The generator $\mathbbm{1}_{4\times 4}$ is, however, not associated with spatial symmetry breaking. We therefore do not take it into account and consider SU($M$) as our starting group instead. that transforms the spinor $\Psi$ from Eq. (6) where, for instance, $M=4$. Moving from sublattice $A$ to $B$ would be described by a matrix $\eta^{12}\in$ SU$(4)$ instead of an element of the Galilei group. The latter is reduced to describe translations with twice the original lattice spacing, $b=2a$. Consequently, translational symmetry breaking on the scale of single plaquettes is now characterized by spontaneous breaking of the subgroup of SU$(4)$ that is compatible with the initial interaction given in Tab. 1.

Let us make this more explicit for the present situation. The compact Lie group SU$(4)$ is generated by a set of 15 generalized Gell-Mann matrices $\{\hat{\lambda}_{1},...,\hat{\lambda}_{15}\}$, see App. B for an explicit representation. The initial interaction, given by the contributions $\sim J_{1,\Lambda}$ in Tab. 1 is invariant under a subgroup, generated by seven of these,

[TABLE]

All of these generators may in principle be subject to spontaneous breaking. Consider now a homogeneous ansatz for the order parameter, i.e. the one corresponding to the BZA phase, in the four-sublattice description,

[TABLE]

This ansatz does not commute with any of the generators given in Eq. (9), separately. However, the combination

[TABLE]

leaves the interaction as well as the order parameter, Eq. (10), invariant. It thus generates an unbroken subgroup.

Indeed, when implementing the special configuration, Eq. (10), corresponding to $Q^{(1)}=Q^{(2)}=Q^{(3)}=Q^{(4)}$, as an initial condition for the RG flow equations, one finds that the divergence is not regularized and solving the flow all the way through $\Lambda\rightarrow 0$ is not possible. The conclusion is that the physics of the $J_{1}$ Heisenberg model requires the subgroup generated by $\lambda_{\pi}$ to be broken as well. To find an appropriate order parameter, all one has to do is to find generators of SU(4) that do not commute with $\lambda_{\pi}$ and add at least one of them to the ansatz from Eq. (10). Eligible candidates are the Gell-Mann matrices

[TABLE]

Not surprisingly, they correspond to modifications of the real or imaginary parts of the initial values of $Q^{(1..4)}$, respectively. Thus, each ansatz with expansion coefficients $Q_{i}\in\mathbbm{R}$

[TABLE]

leads to a maximally broken phase where the RG flow can be expected to converge and (almost) any random ansatz generates a valid flow in the infrared, as observed in Sec. II.2 .

In order to understand the physical meaning of the generators, Eq. (12), recall that the initial interaction, Eq. (1), in the pseudofermion formulation is invariant under an artificial local U(1) symmetry. One finds that the generators $\hat{\lambda}_{i}$ inside of the fermionic order parameter structure

[TABLE]

can be converted into each other by such gauge transformations. The notion of translational symmetry breaking in this context should therefore always be understood modulo the application of such local transformations Affleck et al. (1988).

From the algebraic considerations alone, it is not obvious why the very general ansatz, Eq. (13), should always lead to a $\pi$-flux phase in the infrared. This must indeed be shown by actually solving the flow equations. However, two important conclusions can be drawn already: Firstly, the sublattice ansatz allows to detect finite-ranged spatial symmetry breaking patterns by the familiar divergences. The range can be extended systematically by using larger numbers of sublattices. Secondly, the (non-)occurrence of divergences for specific symmetry breaking or preserving ansätze conveys information about which state is energetically preferred and thus closer to the actual ground state of the system. It thus makes a direct comparison of free energies unnecessary. To emphasize this key feature of our FRG approach, we reiterate that solving the flow equations therefore automatically selects the symmetry-breaking pattern and identifies the energetically favored ground state. The reliability of our approach strongly benefits from this particular finding, as our present truncation scheme is not specifically tailored to compute precise values of the free energy.

III.3 Frustration at large $N$

We now include magnetic frustration and consider the $J_{1}$-$J_{2}$ model for a broad range of the ratio $g=\frac{J_{2}}{J_{1}}$. To that end, we take into account the full initial action and order parameter structure given in Tabs. 1 and 2. The number of couplings generated during the flow increases, too, summing up to 48 differential equations to be solved. The resulting finite-temperature phase diagram is shown in Fig. 4. For sufficiently low temperature, an ordered phase is found for all $g$. While $Q^{(1..4)}$ dominate for $g<1$, finite $Q^{(5,6)}$ become mandatory to describe the phase at $g>1$. The former exhibits $\pi$-flux characteristics, while $Q^{(5,6)}$ do not entail a complex phase for next-nearest neighbor plaquettes. This particular finding is most likely due to our limited spatial range: with four sublattices, no next-to-nearest neighbor $\pi$-flux phase can be described. We would expect to find such a flux phase for $g>1$ as well, provided the number of sublattices would be sufficiently increased. Since the critical temperature does not depend on the spatial structure of the order parameter in this system Affleck and Marston (1988); Arovas and Auerbach (1988), we can still expect to obtain qualitatively correct results.

For $g>1$ and low temperatures, we find a regime of coexistence and a competition of ordering tendencies is indicated, cf. Figs. 4 and 5, by the suppression of the critical temperature or the order parameter in the coexistence region. Indeed, frustration is expected to occur due to a finite $J_{2}$. However, it does not seem sensible to directly relate this finding to the suspected spin liquid regime of the physical model at $N=2$: Firstly, it is widely agreed Chandra and Doucot (1988); Einarsson and Schulz (1995); Schulz et al. (1996); Capriotti and Sorella (2000); Sirker et al. (2006); Mambrini et al. (2006); Richter and Schulenburg (2009); Beach (2009); Jiang et al. (2012); Mezzacapo (2012); Hu et al. (2013); Ren et al. (2014); Gong et al. (2014); Metavitsiadis et al. (2014); Richter et al. (2015); Chou and Chen (2014); Wang et al. (2016); Aktersky and Syromyatnikov (2016); Wang and Sandvik (2018); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Liu et al. (2018) that the latter is situated at about $0.4\lesssim g\lesssim 0.6$ which does not even overlap with the large-$N$ coexistence region. Secondly, it is doubtful whether frustration is an appropriate term to describe the mechanism behind this coexistence. At $N\rightarrow\infty$, substantial contributions of the interactions are suppressed Buessen et al. (2018b), notably those being responsible for magnetic ordering tendencies Baez and Reuther (2017). Those very interactions, commonly referred to as large-$S$ enhanced, are known to facilitate frustration Chandra and Doucot (1988). In our case, instead of being cancelled by (geometrically) competing spin-spin interactions, magnetic ordering is thus suppressed artificially, giving way to a competition between spin-liquid phases instead.

If the mechanism behind the putative spin-liquid phase at $N=2$ is indeed a suppression of magnetic ordering due to frustration and large-$N$ results become a good approximation in this regime, we would therefore expect a nearest-neighbor dominated spin-liquid phase like a $\pi$-flux phase to occur.

IV Phase diagram at $N=2$

At $N=2$, additional fluctuations occur which are absent in the large-$N$ limit. After discussing the implementation of the fermion number constraint, we therefore extend the truncation by taking into account further order parameter structures and newly generated four-fermion interactions. Then, we calculate the phase diagram of the $J_{1}$-$J_{2}$ model for $g\in[0,1]$ with our approach and exhibit the suppression of magnetic correlations upon inclusion of a minimal frequency-dependent self-energy $\Gamma^{(\gamma)}$. In fact, we show how $\Gamma^{(\gamma)}$ leads to a phase diagram which is consistent with the results from previous pf-FRG studies. To improve the understanding of the role of $\Gamma^{(\gamma)}$, we briefly revisit the large-$N$ case and show that, here, the spin-liquid state distinctly originates from the flow of the four-fermion vertex and not from a frequency-dependent self-energy. Including an artificial $\Gamma^{(\gamma)}$ at large $N$ leads to a strong suppression of the spin-liquid state, too. Based on this consideration, we return to $N=2$ and argue that the inclusion of frequency-dependent self-energy induces a damping of both, magnetic as well as bilinearly-ordered spin-liquid states, and therefore does not support the formation of the latter.

IV.1 Fermion number constraint

For large $N$, it is sufficient to implement the fermion number constraint (4) on average Arovas and Auerbach (1988), which greatly simplifies the analysis. Moving to $N=2$, the problem has to be dealt with more carefully. There are at least two viable options on how to faithfully implement the constraint in this situation.

The first and probably better known one is to make use of an artificial local SU(2) symmetry of the model, Eq. (1), introduced due to the pseudofermion formulation Affleck et al. (1988). The constraint as expressed in Eq. (4) explicitly breaks this symmetry down to a local U(1). However, implementing the constraint in the form of a Lagrange multiplier field $\vec{A}$ restores the local SU(2),see Ref. Affleck et al., 1988 for details. Moreover, the local SU(2) now becomes a true gauge symmetry as the Lagrange multiplier plays the role of a gauge field. The confinement physics associated with this gauge field is often taken as a hallmark of spin liquid behavior Savary and Balents (2017). Although the predictive power associated with the gauge field is persuasive, it entails technical challenges that are beyond the scope of this work 222 Dealing with (non-abelian) gauge fields in general is not a easy task within the FRG formalism. It could at best be done in an approximate way and it is not clear in which way a given approximation scheme would influence the fulfillment of the constraint and thus the physics outcome. Furthermore, the field $\vec{A}$ is expected to acquire its own dynamics during the flow. Since we are considering a lattice system, such a gauge field is generally compact and a proper description might have to include vortex or monopole degrees of freedom. Although there is work in this direction in the context of BKT transitions Machado and Dupuis (2010); Krieg and Kopietz (2017), this represents a formidable task on its own and is therefore beyond the scope of this paper.. Instead, we follow a different approach, here.

The other realization of the constraint was put forward by Popov and Fedotov Popov and Fedotov (1988), who showed that equipping the pseudo-fermion Hamiltonian with an imaginary chemical potential

[TABLE]

leads to an exact cancellation of all unphysical contributions in the partition function. The constraint is thus implemented exactly on the level of the microscopic action (5). This procedure is suitable for our FRG approach and the ansatz for the bilinear part of the effective average action $\Gamma_{\Lambda}^{(2\Psi)}$ without order parameters becomes

[TABLE]

with $\mu_{\Lambda_{\rm{UV}}}=\mu_{\mathrm{PF}}$. Here, $\vartheta_{\Lambda}^{-1}$ is the formal inverse of the multiplicative regulator function Enss (2005)

[TABLE]

There is a number of things to be observed about this way of implementing the constraint. Firstly, it is not equivalent to a functional $\delta$ function enforcing Eq. (4) as in the case of the Lagrange multiplier. Hence, it is not permissible to simply replace any $f^{\dagger}_{i\alpha}f^{\phantom{\dagger}}_{i\alpha}$ with the number 1 in the action. Consequently, the structure of the initial interaction is now more complicated than the one given in Eq. (5), because

[TABLE]

and the last term on the right-hand side, corresponding to a density-density interaction, cannot be neglected anymore.

Secondly, the imaginary chemical potential in the ansatz Eq. (16) is equipped with an explicit $\Lambda$ dependence, its initial value set to $\mu_{\mathrm{PF}}$. This may seem odd at first glance, since the constraint (4) is an exact relation that does not seem to lend itself to a continuous RG flow. One should note, however, that the proof of $\mu_{\mathrm{PF}}$ implementing the constraint was given for the microscopic Hamiltonian, i.e. at the scale $\Lambda=\Lambda_{\mathrm{UV}}$. The value of the partition function does, of course, not change under an (exact) RG transformation, but the shape of the Hamiltonian or the action does. Physically, this means that the $f^{(\dagger)}$ fields are systematically dressed. To fulfill the constraint at scales $\Lambda\neq\Lambda_{\mathrm{UV}}$, it is therefore mandatory to let $\mu_{\Lambda}$ evolve with the rest of the effective average action.

Our approach to the solution of the FRG equation (2) is approximate and so is the renormalization of $\mu_{\Lambda}$. Therefore, the implementation of the constraint is not exact for any $\Lambda<\Lambda_{\mathrm{UV}}$, just as for a hypothetical gauge field realization. On the other hand, we do have better control about the approximation error since we keep the truncation exactly at the same level as for the rest of the action. Let us illustrate this point with a toy model of our system. Consider $N=2$ in the limit $J_{2}=0$ without initial order parameters on two sublattices, i.e. $\Psi$ is a two-dimensional spinor. The latter is necessary for an algebraic discrimination between the two contributions to the right hand side of Eq. (18) in momentum space. The four-fermion terms that are present initially or will be generated during the flow are shown in Tab. 3.

The flow equations for the couplings $J_{\Lambda}$, $L_{\Lambda}$, $X_{\Lambda}$ and $\mu_{\Lambda}$ can be found in App. C. Depending on the temperature, the flow is intercepted by a divergence at some scale $\Lambda_{\mathrm{div}}$. While this scale itself does not bear direct physical meaning, its existence signals the onset of ordering Braun (2012). The highest temperature where $\Lambda_{\mathrm{div}}>0$ will be dubbed instability temperature. In Fig. 6, the value of $\Lambda_{\mathrm{div}}$ is plotted for a number of different approximation schemes: no constraint ($\mu_{\Lambda}=0$, dotted line), static chemical potential ($\mu_{\Lambda}=\mu_{\mathrm{PF}}$, dashed line), running chemical potential with (solid line) or without (dot-dashed line) Katanin-feedback of $\mu_{\Lambda}$ on the four-fermion couplings. The impact of these approximations on the instability temperature is quite dramatic and it is therefore vital to include $\mu_{\Lambda}$ into a consistent approximation scheme 333In previous work Roscher et al. (2018) we suggested a scheme that would correspond to the case of a static $\mu_{\Lambda}$. In light of the present findings, we now rather advocate the complete Katanin scheme (solid line in Fig. 6) as the latter is used for the flow with order parameters as well..

One last remark is due to the zero-temperature limit of the problem. By definition in Eq. (15), $\mu_{\mathrm{PF}}$ vanishes for $T\rightarrow 0$. It is not obvious, whether this limit commutes with the path integral, i.e. the computation of the partition function. At least in our toy model, this seems to be the case, though. For all approximations, $\Lambda_{\mathrm{div}}$ smoothly approaches the naïve zero-temperature value.

IV.2 Extended truncation

In the frustrated regime of the $N=2$ case, it is insufficient to consider only the order parameters $Q^{(1..4)}$ and $Q^{(5,6)}$ due to the presence of additional fluctuations away from the large-$N$ limit. Therefore, we have to introduce additional pairing order parameters $\Delta$, allowing us to investigate more classes of possible spin-liquid phases. Our ansatz for the bilinear part $\Gamma_{\Lambda}^{(2\Psi)}$ of the effective average action becomes

[TABLE]

Here, $\mathfrak{Q}(p)$ and $\mathfrak{D}(p)$ are the respective geometric momentum structures. The pairing order parameters are associated with a non-trivial structure $\sim\epsilon_{\alpha\beta}$ in spin space preserving the global SU(2) symmetry. Instead of the totally antisymmetric tensor $\epsilon_{\alpha\beta}$, Pauli matrices $\sigma^{x,z}_{\alpha\beta}$ or linear combinations would be permissible as well. Due to SU(2) symmetry, however, this does not change the structure of the flow equations.

Introducing the new order parameters has consequences for the four-fermion sector of the model as well. For $\Delta^{(j)}=0$, Eq. (7) is the only possible interaction structure, different combinations of $\eta^{X/Y}$ determining the number of couplings to be taken into account. Once any $\Delta^{(j)}>0$, there are five new structures that can and generally will be generated during the flow, see Tab. 4. Note that, SU(2) symmetry keeps the spin structure of the interactions simple in the sense that $\delta_{\alpha\beta}$ and $\epsilon_{\alpha\beta}$ are associated to fixed Nambu sectors during the entire RG flow.

As opposed to the large-$N$ case, all combinations of $\eta^{X/Y}$ matrices can now contribute to the RG flow. Taking into account all of these combinations for all Grassmann structures would correspond to a total of 1536 couplings in the four-fermion sector. There is, however, a number of symmetries that can be exploited: for instance, the two mixed Grassmann structures are hermitean conjugates of each other and thus do not contain independent information. Furthermore, the choice of the $\eta^{X/Y}$ matrices themselves is constrained by hermiticity as well. Taking all of this into account reduces the number of independent (complex) couplings from the four-fermion sector down to 624. The flow equations for these are then solved by step-size-controlled numerical integration.

When translating the initial action (5) into momentum space, the definition of nearest or next-nearest neighbor coupling guarantees a unique representation of the vertex’ momentum structure with respect to the choice of sublattices. In contrast to the large-$N$ case, this structure is not invariant under the RG flow at $N=2$ anymore. Longer-ranged interactions are generated. Whenever the range of such an interaction exceeds the enlarged unit cell given by the sublattice representation, the momentum structure is altered. In order to properly account for this effect, projection onto the momentum structure itself would be required. Here, we strictly truncate our interaction vertices to take into account the nearest-possible neighbor sites, only. In practice, this corresponds to the following generic momentum structure

[TABLE]

The $k_{i}$ are the momenta of the vertex’ fields. $\{a_{x/y},b_{x/y},c_{x/y}\}$ are determined by $\eta^{X/Y}$, i.e. the geometry of the vertex and $\{i,j,k\}$ vary with the Grassmann sector considered, see Tab. 4. For $N\rightarrow\infty$, we recover the momentum structures employed in Sec. II.2.

We do not expect this truncation to interfere with the regularization of the RG flow through emergent spin-liquid order. This is in analogy to the analysis presented in Sec. III.2: there, we showed that, in the two-sublattice parametrization, the uniform BZA state regularizes the flow, but in the four-sublattice parametrization the additional translational symmetry breaking of the $\pi$-flux state is exhibited. The level of truncation was equivalent to eq. (20) in both instances. We expect this mechanism to extend to parametrizations with more sublattices.

IV.3 Phase diagram & frequency dependence

We first discuss our numerical findings for $J_{2}=0$. In this case, a Néel-ordered magnetic phase is expected to occur at sufficiently low temperatures Manousakis (1991). To compute infrared properties of such a phase in our approach, one would have to introduce magnetic order parameters, i.e. terms $\sim M^{(n)}\Psi^{\dagger}_{\alpha}\tau_{\alpha\beta}\eta^{X}\Psi_{\beta}$ where $\tau_{\alpha\beta}\neq\delta_{\alpha\beta}$ in spin space. Here, we are not interested in magnetic orders and we refrain from taking them into account. This greatly reduces the numerical cost. For $T/J_{1}<0.82$ we find a divergence of the four-fermion couplings that does not trigger a commensurate growth of any of the viable bilinear spin-singlet orders $Q^{(i)}_{\Lambda}$ or $\Delta^{(j)}_{\Lambda}$. Therefore, the divergence indicates a spontaneous breaking of the global SU(2) and hence magnetic ordering.

Considering $J_{2}>0$, the phase diagram in Fig. (7) ensues (solid line), where magnetic ordering occurs for all $g=J_{2}/J_{1}$ with just a little reduction of the critical temperature around the suspected regime of frustration, $g\sim 0.5$. This gross overestimation of magnetic ordering tendencies is expected in the presented static approximation, where explicit frequency dependences of the four-fermion vertices or self-energy contributions have been neglected: properly taking frequency dependence into account within the pf-FRG scheme Reuther and Wölfle (2010) leads to a magnetic phase diagram which is well compatible with the bulk of the literature, in particular with regard to the existence of a non-magnetized regime for $0.4\lesssim g\lesssim 0.6$.

Here, we shed light on the mechanism behind the suppression of magnetic ordering aiming at a better characterization of the spin liquid behavior that is facilitated instead: in Ref. Reuther and Wölfle, 2010 it was shown that the phenomenological consequences of frequency dependence, such as finite pseudofermion lifetime, can be reproduced within the static approximation if an ad-hoc self energy

[TABLE]

is added to the ansatz for the effective average action. Being associated to an odd frequency structure, $\gamma$ is not renormalized as long as the vertex functions remain frequency-independent Reuther and Wölfle (2010). By choosing larger values for $\gamma$, magnetic order is more and more suppressed until, eventually, the regime around $g\sim 0.5$ is indeed free of magnetic instabilities all the way down to $T=0$, cf. Fig. 7.

IV.4 Large-$N$ spin liquid & frequency dependence

To assess the influence of frequency dependence on possible spin liquid phases, let us return for a moment to the large-$N$ situation: There, no change of the critical temperature is found when comparing a fully frequency dependent pf-FRG scheme to the simplified one at hand Buessen et al. (2018b); Roscher et al. (2018). The reason for this behavior can be understood from perturbative diagrammatics. An imaginary self-energy like (21) can naturally not occur from a one-loop diagram as the initial value for the four-fermion vertex is frequency independent. At two loop order, the diagram given in Fig. 8a accounts for a non-trivial renormalization of the frequency dependence of the self-energy. In our FRG flow, it is effectively generated by a dressed vertex that is fed into the flow of the self energy. Since only one of the two loops contributes a factor of $N$ due to a freely summed spin index, the overall value of the diagram is $\sim\frac{1}{N}$ and thus suppressed for $N\rightarrow\infty$. Thus, even if the frequency dependence of vertices was taken into account, no such dependence of the self energy can be generated at large $N$. It is therefore safe to say that in the present system, the occurrence of spin-liquid phases at large $N$ is solely due to the peculiar structure of the four-fermion vertex flow: only the diagram with a particle-hole loop (Fig. 8b) contributes while, in particular, the RPA-type diagram (Fig. 8c), which is typically associated with magnetic ordering, is suppressed.

IV.5 $N=2$ case

Considering the case $N=2$, both of these two-particle diagram types and the self energy (21) become important. We do know the effect of finite damping $\gamma>0$ on magnetic ordering (see Fig. 7), but its impact on potential spin liquid phases is still unclear, yet of utmost importance for the understanding of the highly frustrated $g\sim 0.5$ regime. To shed light on this question, we need to investigate the interplay of $\Gamma_{\Lambda}^{(\gamma)}$ and spin-liquid favoring fluctuations represented by the diagram in Fig. 8b. To this end, we artificially introduce the self energy (21) into the large-$N$ flow equations of the $J_{1}$ model for two sublattices, which are by construction dominated by this type of fluctuations. This gives us a handle on the influence of $\gamma>0$ on spin liquid ordering tendencies also beyond the large-$N$ limit. The so-obtained results are shown in Fig. 9.

We find that bilinear spin liquid ordering is suppressed upon introduction of an artificial damping $\gamma$ as well. In fact, at a value of $\gamma=0.3$, no such phase could be discerned anymore. This is particularly noteworthy, since for $\gamma=0.3$, the suspected spin liquid regime around $g\sim 0.5$ at $N=2$ (cf. Fig. 7, dotted line) is still occupied by magnetic phases. Running the complete $N=2$ flow equations again results in a divergence, no matter if the order parameters from Eq. (19) are present or not.

We draw two conclusions from these results. Firstly, frequency dependence alters quantitative features such as the values of critical temperatures or the position of phase boundaries. However, since damping as its main consequence affects both, conventional (magnetic) as well as bilinear spin liquid ordering, we do not expect the appearance of physical spin liquid phases which are not present in the static FRG approximation. This solidifies the reliability of our approach. Secondly, the primary goal of this work was to shed some further light on the nature of the frustrated regime around $g\sim 0.5$. Importantly, we can conclude that our approach does not support the emergence of any particular type of bilinearly ordered spin liquid, since we included all possible bilinear order parameters that do not break the global SU(2) symmetry. This indicates that bilinear spin-liquid order is not an appropriate description of the physical phenomena occurring in the $J_{1}$-$J_{2}$ model in the regime $0.4\lesssim g\lesssim 0.6$.

By construction, our approach is sensitive to the occurrence of any bilinear order parameter, signaled by a divergence in the four-fermion sector of the effective average action. However, ordering tendencies that need to be described by expectation values of more fermion fields are not seen at the given truncation level. The results presented, here, therefore strongly hint towards such a more complicated structure of the spin liquid phase. Possible candidates are dimer valence bond solids Sirker et al. (2006); Metavitsiadis et al. (2014); Haghshenas and Sheng (2018) (involving four fermion fields) or a four-spin (eight-fermion) plaquette order Capriotti and Sorella (2000); Mambrini et al. (2006); Gong et al. (2014); Aktersky and Syromyatnikov (2016).

V Conclusions & Outlook

To summarize, we have introduced a cluster pf-FRG approach and applied it to the $J_{1}$-$J_{2}$ Heisenberg model on the square lattice. First, a large-$N$ analysis indicated the necessity to incorporate spatially structured order parameters by means of the cluster extension. Moving to the physical case at $N=2$, we carefully analyzed the implementation of the fermion number constraint using an imaginary chemical potential according to the Popov-Fedotov method. We added a minimal frequency-dependent self-energy to our truncation to incorporate the effects of finite pseudofermion lifetime. By doing so we were able to reproduce the non-magnetic phase at intermediate coupling ratio $0.4\lesssim g\lesssim 0.6$. However, no signs of bilinear spin-liquid order emerged in this regime. In fact, any instability encountered in this model could reliably be related to magnetic ordering tendencies instead. This suggests that the low temperature regime of the strongly frustrated $J_{1}$-$J_{2}$ Heisenberg model is occupied by a state that needs to be characterized by the expectation value of at least four fermion fields, the simplest incarnation of which includes valence bond solid or plaquette ordered states, but also more complex spin liquid states.

There are multiple directions for future work. As indicated in Sec. II.2, the sublattice representation opens up our method for further, more complicated lattice geometries. A natural next step is to perform an analogous analysis for the Kitaev model on the honeycomb lattice. While the exact solution is formulated in terms of Majorana fermions Kitaev (2006), it has been shown that an equivalent pseudofermion representation exists Burnell and Nayak (2011). In this case, the QSL phase is indeed characterized by the order parameters introduced in Eq. (19). Having reproduced these findings, one should be able to analyze extensions of the Kitaev model with only minor additional effort.

Acknowledgements.

We thank F. L. Buessen and D. Kiese for ongoing collaborations on related projects and comments on the manuscript. We acknowledge valuable discussions with J. Reuther and B. Sbierski and their comments on the manuscript. We also thank J. Braun and N. Dupuis for insightful discussions. This work was partially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 277146847 – SFB 1238 (project C03). S.D. and D.R. acknowledge support by the German Research Foundation (DFG) through the Institutional Strategy of the University of Cologne within the German Excellence Initiative (ZUK 81).

Appendix A Large-$N$ structure of the fluctuation matrix

The right hand side of the flow equation (2) is determined by the second functional derivative of the effective average action, $\Gamma_{\Lambda}^{(2)}=\mathcal{P}_{\Lambda}+\mathcal{F}_{\Lambda}.$ Here, $\mathcal{F}_{\Lambda}$ is the field-dependent part of $\Gamma_{\Lambda}^{(2)}$, called fluctuation matrix. Applying the point-like projection $\Psi_{\alpha,\omega_{m}-\nu_{n},\mathbf{p-q}}=\Psi_{\alpha}\delta_{m,n}\delta(\mathbf{p-q})$, its value with respect to the generic interaction (7) is given by

[TABLE]

where

[TABLE]

The flow equation (2) may be expanded in powers of $\mathcal{F}_{\Lambda}$. The quadratic contribution

[TABLE]

of this expansion determines the flow of the fermionic interaction terms themselves. Thus, any $\partial_{\Lambda}J_{n}\sim\frac{1}{N}$. Only terms where the trace operation on the right hand side of the Wetterich equation contributes an extra factor of $N$ avoid being suppressed in the limit $N\rightarrow\infty$. By construction, $\mathcal{P}_{\Lambda}$ is diagonal in the spin indices. To achieve an overall $N=\mathrm{Tr}[\delta_{\gamma\delta}]$, only components $\sim\delta_{\gamma\delta}$ in eqns. (23) can contribute. These are given in the second lines of eqns. (23b) and (23c). Multiplying those terms according to eq. (24) can never result in any product of matrices $\sigma^{x}\cdot\eta^{Y}$ sandwiched between $\Psi^{\dagger},\Psi$. Since this is the only way to generate new matrix structures besides the existing $\eta^{X/Y}$, it can be concluded that no such new structures can be generated at large $N$.

Appendix B Gell-Mann matrices of SU(4)

For definiteness, we here list the generators of SU(4) in the fundamental representation as used in Sec. III.2.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Flow equations of the toy model

Employing the regulator (17) the flow equations for the toy model of Sec. IV.1 are given by

[TABLE]

where

[TABLE]

and

[TABLE]

In the limit $T\rightarrow 0$, $\mu_{\Lambda}$ vanishes. The flow equations thus simplify dramatically, yielding

[TABLE]

We use these equations to compute the $T=0$ data point in Fig. 6.