Classical spiral spin liquids as a possible route to quantum spin liquids

TL;DR

This paper explores how classical spiral spin liquids on bipartite lattices can evolve into quantum spin liquids when quantum fluctuations are considered, using theoretical methods to identify potential quantum phases.

Contribution

It introduces a general scheme to construct classical spiral spin liquids on bipartite lattices and investigates their quantum counterparts using PFFRG, revealing possible routes to quantum spin liquids.

Findings

Classical spiral spin liquids can evolve into quantum spin liquids with quantum fluctuations.

Remnants of spiral degeneracies persist in quantum $S=1/2$ regimes.

The study provides a framework for designing new quantum spin-liquid candidates.

Abstract

Quantum spin liquids are long-range entangled phases whose magnetic correlations are determined by strong quantum fluctuations. While an overarching principle specifying the precise microscopic coupling scenarios for which quantum spin-liquid behavior arises is unknown, it is well-established that they are preferably found in spin systems where the corresponding classical limit of spin magnitudes $S\rightarrow\infty$ exhibits a macroscopic ground state degeneracy, so-called classical spin liquids. Spiral spin liquids represent a special family of classical spin liquids where degenerate manifolds of spin spirals form closed contours or surfaces in momentum space. Here, we investigate the potential of spiral spin liquids to evoke quantum spin-liquid behavior when the spin magnitude is tuned from the classical $S\rightarrow\infty$ limit to the quantum $S=1/2$ case. To this end, we first use the Luttinger-Tisza method to formulate a general scheme which allows one to construct new spiral spin liquids based on bipartite lattices. We apply this approach to the two-dimensional square lattice and the three-dimensional hcp lattice to design classical spiral spin-liquid phases which have not been previously studied. By employing the pseudofermion functional renormalization group (PFFRG) technique we investigate the effects of quantum fluctuations when the classical spins are replaced by quantum $S=1/2$ spins. We indeed find that extended spiral spin-liquid regimes change into paramagnetic quantum phases possibly realizing quantum spin liquids. Remnants of the degenerate spiral surfaces are still discernible in the momentum-resolved susceptibility, even in the quantum $S=1/2$ case. In total, this corroborates the potential of classical spiral spin liquids to induce more complex non-magnetic quantum phases.

Equations

$$H=\frac{1}{2}\sum _{ij}J({\mathbf{R}}_{ij})\mathbf{S}({\mathbf{r}}_i)\mathbf{S}({\mathbf{r}}_j).$$

$${\tilde{\mathbf{S}}}_{\alpha }(\mathbf{q})$$

$${\tilde{J}}_{\alpha \beta }(\mathbf{q})$$

$$H=\sum _{\mathbf{q}}\sum _{\alpha ,\beta }{\tilde{J}}_{\alpha \beta }(\mathbf{q}){\tilde{\mathbf{S}}}_{\alpha }(\mathbf{q}){\tilde{\mathbf{S}}}_{\beta }(-\mathbf{q}),$$

$${\tilde{\mathbf{S}}}_{\alpha }(\mathbf{q})=\sum _{\nu =1}^{\Gamma }{\mathbf{w}}_{\nu }(\mathbf{q})u_{\alpha }^{\nu }(\mathbf{q})$$

$$H=\sum _{\mathbf{q}}\sum _{\nu =1}^{\Gamma }{\lambda }_{\nu }(\mathbf{q}){\left|{\mathbf{w}}_{\nu }(\mathbf{q})\right|}^2.$$

$$|{\mathbf{S}}_i{|}^2=1\forall i,$$

$$\sum _{i=1}^N|{\mathbf{S}}_i{|}^2=N.$$

$$\mathbf{S}({\mathbf{r}}_{i\in \alpha })=\frac{1}{\sqrt{N/\Gamma }}\sum _{\mathbf{q}}{e}^{i{\mathbf{qr}}_i}{\tilde{\mathbf{S}}}_{\alpha }(\mathbf{q}),$$

$$\sum _{\mathbf{q}}\sum _{\nu =1}^{\Gamma }{\left|{\mathbf{w}}_{\nu }(\mathbf{q})\right|}^2=N.$$

$$H\ge N{\lambda }_{LT},$$

$${\left|{\mathbf{w}}_{\nu }(\mathbf{q})\right|}^2=\{\begin{array}{ll}\frac{N}{2}& \text{if }\mathbf{q}=\pm {\mathbf{q}}_{LT}\wedge \nu ={\nu }_{LT}\\ 0& \text{otherwise}\end{array}.$$

$$\mathbf{S}({\mathbf{r}}_{i\in \alpha })=\frac{1}{\sqrt{2\Gamma }}\sum _{\mathbf{q}=\pm {\mathbf{q}}_{\mathbf{LT}}}\hat{\mathbf{n}}(\mathbf{q})u_{\alpha }^{LT}(\mathbf{q})e^{i{\mathbf{qr}}_i}$$

$$J_{\alpha \alpha }(\mathbf{q})=J_2^{\alpha }\sum _{i=1}^n\cos (\mathbf{q}{\mathbf{a}}_i)\equiv J_2^{\alpha }f(\mathbf{q}).$$

$$J_{\alpha \beta }(\mathbf{q})=\frac{J_1}{2}\sum _{{\mathbf{R}}_{\alpha \beta }^{\text{nn}}}{e}^{i\mathbf{q}{\mathbf{R}}_{\alpha \beta }^{\text{nn}}},$$

$$\underline{\mathbf{J}}(\mathbf{q})=\left(\begin{array}{cc}{J}_{11}(\mathbf{q})& J_{12}(\mathbf{q})\\ J_{12}^{\ast }(\mathbf{q})& J_{11}(\mathbf{q})\end{array}\right).$$

$${\mathbf{u}}_{\nu }(\mathbf{q})$$

$${\lambda }_{\nu }(\mathbf{q})$$

$$e^{i\varphi (\mathbf{q})}=\sqrt{\frac{J_{12}(\mathbf{q})}{J_{12}^{\ast }(\mathbf{q})}}.$$

$$\mathbf{S}({\mathbf{r}}_{i\in \alpha })=\left(\begin{array}{c}\cos ({\mathbf{q}}_{LT}{\mathbf{r}}_i+{\varphi }_{\alpha })\\ \sin ({\mathbf{q}}_{LT}{\mathbf{r}}_i+{\varphi }_{\alpha })\\ 0\end{array}\right),$$

$$\nabla {\lambda }_{\nu }(\mathbf{q})=\nabla J_{11}(\mathbf{q})-\frac{1}{2|J_{12}(\mathbf{q})|}\nabla |J_{12}(\mathbf{q}){|}^2=0.$$

$$|J_{12}(\mathbf{q})|=\frac{|J_1|}{2}\sqrt{\left(\sum _{{\mathbf{R}}_{12}^{\text{nn}}}{e}^{i\mathbf{q}{\mathbf{R}}_{12}^{\text{nn}}}\right){\left(\sum _{{\mathbf{R}}_{12}^{\text{nn}}}{e}^{i\mathbf{q}{\mathbf{R}}_{12}^{\text{nn}}}\right)}^{\ast }}.$$

$$|J_{12}(\mathbf{q})|=\frac{|J_1|}{2}\sqrt{z+pg(\mathbf{q})},$$

$$J_2^{\star }\nabla f(\mathbf{q})-\frac{|J_1|}{4}\frac{p}{\sqrt{z+pg(\mathbf{q})}}\nabla g(\mathbf{q})=0.$$

$$J_2^{\star }-\frac{|J_1|}{4}\frac{p}{\sqrt{z+pf(\mathbf{q})}}=0,$$

$$f(\mathbf{q})=p{\left(\frac{J_1}{4J_2^{\star }}\right)}^2-\frac{z}{p}.$$

$$J_{12}(\mathbf{q})=\frac{J_1}{2}{e}^{-i\mathbf{q}\mathbf{h}}\left(1+e^{i\mathbf{q}\mathbf{a}}+e^{i\mathbf{q}\mathbf{b}}+e^{i\mathbf{q}(\mathbf{a}+\mathbf{b})}\right),$$

$$|J_{12}(\mathbf{q})|$$

$$+2\cos (q_a)+2\cos (q_b)]{\}}^{\frac{1}{2}}$$

$$\equiv \frac{J_1}{2}\sqrt{4+2g(\mathbf{q})}.$$

$$f(\mathbf{q})=\frac{1}{8}{\left(\frac{J_1}{J_2^{\star }}\right)}^2-2.$$

$$J_{12}(\mathbf{q})$$

Full text

I Introduction

Quantum spin liquids are prime examples for exotic phases of matter where quantum phenomena and many-body effects combine to create novel emergent properties such as topological order, long-range entanglement and fractional quasiparticle excitations.Anderson (1973); Balents (2010); Savary and Balents (2017) Broadly speaking, quantum spin-liquid behavior may arise in spin systems with sufficiently small magnitudes of the local magnetic moments (typically spin-1/2) when strongly frustrating interactions hinder the system from developing long-range magnetic order in the ground state.

Despite the enormous interest of the condensed matter community in these phases and various decades of intense research activities, quantum spin liquids, however, remain elusive and have so far not been unambiguously detected in real materials. The serious challenges associated with the research on quantum spin liquids also concerns their theoretical understanding; for example, a complete theory describing their correlations, excitations and topological properties does currently not exist. Similarly, a rigorous criterion specifying in which precise microscopic spin models quantum spin liquids exist is also not known. There is, nevertheless, a widespread perception that quantum spin liquids primarily occur in spin systems where the corresponding classical model exhibits a macroscopic ground state degeneracy. A large degeneracy enables stronger thermal and/or quantum fluctuations such that long-range magnetic order may even be suppressed in the classical limit which is known under the term ‘classical spin liquid’.Balents (2010); Lacroix et al. (2011)

Previously, the approach of stabilizing quantum spin liquids starting from classical spin liquids has mostly been pursued for kagomeSachdev (1992); Lecheminant et al. (1997); Yan et al. (2011); Han et al. (2012), pyrochloreRamirez et al. (1999); Bramwell and Gingras (2001); Castelnovo et al. (2008); Gingras and McClarty (2014) and related lattices where corner-sharing triangles or tetrahedra give rise to an ‘ice rule’ constraint which is the origin of a residual entropy at zero temperature. A second and less explored possibility of achieving a macroscopic classical ground state degeneracy is via spiral spin liquids.Bergman et al. (2007); Gao et al. (2016); Iqbal et al. (2018); Buessen et al. (2018a); Mulder et al. (2010); Henley (1987); Ghosh et al. (2019a); Attig and Trebst (2017) In these phases, a special interplay of lattice geometries and frustrating long-range interactions induces an exact degeneracy of spiral ground states which is usually of subextensive type, i.e., in two-dimensional systems the wave vectors defining these spirals form closed contours in momentum space while in three dimensions they form surfaces.

Two paradigmatic systems harboring spiral spin liquids are the Heisenberg models on the two-dimensional honeycomb and the three-dimensional diamond lattice with nearest neighbor $J_{1}$ and antiferromagnetic second neighbor $J_{2}>0$ couplings.Bergman et al. (2007); Gao et al. (2016); Iqbal et al. (2018); Buessen et al. (2018a); Mulder et al. (2010) Even though these systems are bipartite and, hence, show unfrustrated ferromagnetic or Néel magnetic order in the $J_{1}$-only case, when $J_{2}>0$ is tuned beyond a certain critical coupling ratio $J_{2}/|J_{1}|$, spiral contours or surfaces begin to form marking the onset of a spiral spin liquid. Various numerical studies indeed indicate that both models retain their non-magnetic ground states when replacing the classical $S\rightarrow\infty$ spins by quantum $S=1/2$ spins possibly realizing a quantum spin-liquid phase.Buessen et al. (2018a); Reuther et al. (2011a); Zhu et al. (2013); Gong et al. (2013); Zhang and Lamas (2013); Li et al. (2012) Despite the fragility of spiral spin liquids with respect to order-by-disorder effects and perturbing longer-range couplings which may easily lift the degeneracy,Bergman et al. (2007); Iqbal et al. (2018); Buessen et al. (2018a) an approximate version of this phase has been experimentally identified in the spin-5/2 diamond lattice compound MnSc2S4.Gao et al. (2016) More recently, experimental and theoretical investigations of the spin-1 diamond spinel NiRh2O4 indicate that this material might realize a situation where quantum spin-liquid behavior originates from a spiral spin liquid in the classical $S\rightarrow\infty$ limit.Buessen et al. (2018a); Chamorro et al. (2018) One reason why this scenario of stabilizing a quantum spin liquid is less explored as compared to the aforementioned kagome and pyrochlore geometries is because there is currently no general criterion known which allows one to systematically construct models with spiral degeneracies.

The purpose of this work is two-fold. Firstly, based on the classical Luttinger-Tisza method,Luttinger and Tisza (1946); Luttinger (1951) we formulate an approach to create new models for spiral spin liquids on bipartite lattices. The known spiral spin liquids on the honeycombMulder et al. (2010), diamondBergman et al. (2007) and bcc latticesAttig and Trebst (2017) all fall into the category of systems that can be described within this approach. We then extend these considerations and construct two more models harboring spiral spin liquids which are less explored or have not been studied before. The first is a Heisenberg model on the two-dimensional square lattice with antiferromagnetic nearest neighbor $J_{1}$, second neighbor $J_{2}$ and third neighbor $J_{3}$ interactions. Interestingly, even though this system has been extensively studied in the existing literature,Danu et al. (2016); Iqbal et al. (2016a); Chubukov (1991); Ferrer (1993); Mambrini et al. (2006); Arlego and Brenig (2008); Reuther et al. (2011b); Ralko et al. (2009) the spiral spin liquid occurring for $J_{2}/J_{1}>1/4$ and $J_{3}/J_{2}=1/2$ has so far rarely been discussed.Ceccatto et al. (1993); Moreo et al. (1990) Additionally, we construct a spiral spin liquid on the three-dimensional hcp lattice (hcp stands for ‘hexagonal close packing’) which requires a total of four antiferromagnetic couplings. We show that there is a one-parameter manifold of couplings (i.e. two coupling ratios are fixed) for which the system develops degenerate spiral surfaces in momentum space.

Our second objective is to demonstrate that these spiral spin liquids are good candidate systems for realizing a quantum spin liquid in the $S=1/2$ case. Therefore, we employ the pseudofermion functional renormalization group (PFFRG) methodReuther and Wölfle (2010); Baez and Reuther (2017) which is capable of treating strongly frustrated spin systems even in the case of complex two-dimensionalReuther et al. (2011a); Reuther and Wölfle (2010); Reuther et al. (2011b); Baez and Reuther (2017); Iqbal et al. (2016a); Balz et al. (2016); Keleş and Zhao (2018) or three-dimensionalIqbal et al. (2019, 2018); Buessen et al. (2018a); Iqbal et al. (2017, 2016b) lattice geometries and in the presence of long-range couplings.Reuther et al. (2011a, b); Keleş and Zhao (2018); Iqbal et al. (2016a, 2018); Buessen et al. (2018a); Iqbal et al. (2017, 2016b) Our numerical results indicate that for both models the coupling regimes of classical spiral degeneracies indeed host extended non-magnetic phases in the $S=1/2$ case. The weight distributions in the static but momentum-resolved spin susceptibilities also show that residual features of the degenerate spiral surfaces still survive in the quantum case, however, depending on the precise coupling parameters the sizes of these surfaces are seen to deviate from their classical values. In total, the two models which we investigate in the classical and in the quantum case constitute additional examples demonstrating that classical spiral spin liquids provide a promising platform for the search for quantum spin liquids.

The remainder of this work is organized as follows: In Sec. II we investigate under which conditions spiral degeneracies occur on bipartite lattices. These considerations are based on the classical Luttinger-Tisza approach which is briefly introduced in Sec. II.1. The general scheme for designing spiral spin liquids on bipartite lattices is developed and formulated in Sec. II.2. The following Sec. III makes this procedure more concrete by demonstrating its capability for the square (Sec. III.1) and for the hcp lattice (Sec. III.2). The effects of quantum fluctuations on both systems are studied in Secs. IV.1 and IV.2, respectively, where the PFFRG method is applied to calculate momentum-resolved susceptibilities and to map out phase diagrams. We conclude the paper in Sec. V with a short summary of the main results and an outlook for future directions of research.

II Ground state degeneracies in classical bipartite Heisenberg systems

In this section we apply the classical Luttinger-Tisza method Luttinger and Tisza (1946); Luttinger (1951) to find a criterion for the formation of spiral spin liquids that can be used for the identification of new systems with this property. Since the Luttinger-Tisza method will be essential for formulating this criterion, we will first give a brief introduction into this approach. The considerations below crucially rely on the fact that our models are defined on bipartite lattices.

II.1 The Luttinger-Tisza method

In the classical Heisenberg model, spins are represented as three-dimensional vectors, such that the Hamiltonian of a system with $N$ classical spins $\bm{S}_{i}=\bm{S}(\bm{r}_{i})$ coupled by exchange interactions $J_{ij}=J(\bm{R}_{ij})$ can be written as

[TABLE]

Here, $\bm{r}_{i}$ is the real space position of site $i$ and $\bm{R}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ is the distance between two sites $i$ and $j$. Even in the classical case, finding the ground state of Eq. (1) is usually a non-trivial minimization problem. The Luttinger-Tisza method solves this problem, at least approximatively, but can also be exact in special cases, as will be discussed in more detail below. This approach starts by defining the Fourier transforms of the spins $\tilde{\bm{S}}_{\alpha}(\bm{q})$ and the exchange interactions $\tilde{J}_{\alpha\beta}(\bm{q})$ on each of the sublattices $\alpha$,

[TABLE]

Here, we a have assumed that the lattice has $\Gamma$ sites per unit cell, i.e., $\alpha=1,\ldots,\Gamma$ runs over all the $\Gamma$ sublattices and $\frac{N}{\Gamma}$ is the number of sites in each sublattice. Furthermore, in the definition of $\tilde{J}_{\alpha\beta}(\bm{q})$ the index $i\in\alpha$ denotes an arbitrary but fixed site on sublattice $\alpha$. Note that $\tilde{J}_{\alpha\beta}(\bm{q})$ can be interpreted as the elements of a $\Gamma\times\Gamma$ matrix $\underline{\bm{J}}(\bm{q})$ which contains all the interactions between the sublattices $\alpha$ and $\beta$. Rewriting the Hamiltonian in Eq. (1) in terms of $\tilde{\bm{S}}_{\alpha}(\bm{q})$ and $\tilde{J}_{\alpha\beta}(\bm{q})$ leads to

[TABLE]

where the sum $\sum_{\bm{q}}$ runs over all wave vectors in the first Brillouin zone. It is crucial to realize that the matrix $\underline{\bm{J}}(\bm{q})$ is hermitian and, hence, has real eigenvalues. Its normalized eigenvectors form an orthonormal basis of the $\Gamma$-dimensional vector space. Each cartesian component of the Fourier-transformed spin vectors can thus be expressed as a linear combination of eigenvectors $\bm{u}_{\nu}(\bm{q})$ of $\underline{\bm{J}}(\bm{q})$ and, therefore, we may write

[TABLE]

with $u_{\alpha}^{\nu}(\bm{q})$ being the $\alpha$-th component of the $\nu$-th eigenvector and $\bm{w}_{\nu}(\bm{q})$ the $\nu$-th vector that determines the different cartesian components of $\tilde{\bm{S}}_{\alpha}(\bm{q})$. Inserting this into Eq. (4) results in a Hamiltonian that depends only on the eigenvalues $\lambda_{\nu}(\bm{q})$ and the coefficients $\bm{w}_{\nu}(\bm{q})$,

[TABLE]

Additionally, we have the condition that all spins are normalized,

[TABLE]

which is known as the ‘strong constraint’. Minimizing the Fourier-transformed Hamiltonian in Eq. (6) under the constraint in Eq. (7) does not yet simplify the problem. The key conceptual step proposed by Luttinger and Tisza to approximately solve the minimization problem amounts to replacing the ‘strong constraint’ by the so-called ‘weak constraint’. The latter is much less restrictive, as it only constrains the total spin, instead of imposing separate conditions for all $N$ particles,

[TABLE]

Evidently, all solutions to the problem that satisfy the strong constraint fulfill the weak one as well, whereas the opposite is generally not the case. Therefore, only those solutions found under the weak constraint that additionally satisfy the strong constraint describe physical ground states. In Sec. II.2 it will be shown that under rather mild additional assumptions the Luttinger-Tisza method becomes exact for all bipartite lattices.

Using the inverse Fourier transform of the spins on the $\alpha$-th sublattice,

[TABLE]

the weak constraint can be re-expressed in Fourier space,

[TABLE]

Using Eq. (6) and Eq. (10), a lower limit for the energy follows immediately,

[TABLE]

where the Luttinger-Tisza eigenvalue $\lambda_{LT}\equiv\min{\left\{\lambda_{\nu}(\bm{q})\right\}}$ is defined as the minimum out of the set of eigenvalues of $\underline{\bm{J}}(\bm{q})$ with respect to all wave vectors. The energy can therefore be minimized by setting the coefficients $\bm{w}_{\nu}(\bm{q})$ to be only non-zero for the values $\pm\bm{q}_{LT}$ and $\nu_{LT}$ that minimize $\lambda_{\nu}(\bm{q})$,

[TABLE]

This can be seen by first noting that from Eq. (3) it follows $J_{\alpha\beta}(-\bm{q})=J_{\alpha\beta}^{*}(\bm{q})=J_{\beta\alpha}(\bm{q})$ and thus $\lambda(\bm{q})=\lambda(-\bm{q})$. In the case of several equivalent minima at wave vectors $\{\bm{q}_{LT}\}$, one may chose non-zero $\bm{w}_{\nu}(\bm{q})$ for any $\bm{q}\in\{\pm\bm{q}_{LT}\}$, if the resulting state satisfies the strong constraint. These coefficients by construction satisfy the weak constraint and lead to a spin configuration that is obtained by transforming Eq. (5) back to real space,

[TABLE]

Here, $\hat{\bm{n}}(\bm{q})$ is the unit vector along $\bm{w}_{\nu_{LT}}(\bm{q})$ and $u_{\alpha}^{LT}$ is the $\alpha$-th component of the eigenvector of $\underline{\bm{J}}(\bm{q}_{LT})$ that corresponds to the smallest eigenvalue.

Finally, the physical solutions are only those that satisfy the strong constraint in Eq. (7). It is now obvious why both $+\bm{q}$ and $-\bm{q}$ are taken into account: Since the normalized vector $\hat{\bm{n}}(\bm{q})$ is still undetermined, different complex phases for its $x$- and $y$-components must be chosen such that the spins $\bm{S}_{i}$ are real, which produces coplanar spirals of the form $\bm{S}(\bm{r}_{i})=(\cos{(\bm{q}_{LT}\bm{r}_{i})},\sin{(\bm{q}_{LT}\bm{r}_{i})},0)$. This solution is most evident for mono-atomic Bravais lattices where $u^{LT}(\bm{q})=1$ and the method becomes exact (at least in the spin-isotropic Heisenberg case). For non-Bravais lattices, the ground state cannot generally be represented by such a spiral state. It should further be stressed that any rotation of this spin configuration is also a possible solution as the isotropic Heisenberg model is invariant under uniform rotations of all spins.

If the Luttinger-Tisza method results in one or more possible spin configurations that fulfill the strong constraint, the physical ground state of the system is found. Vice versa, if no solutions can be obtained that have a normalized spin configuration, then the method fails, resulting only in a lower limit for the system’s energy given by Eq. (11). A generalized version of the weak constraint that allows for distinct magnitudes of spins on each sublattice has been proposed by D. H. Lyons and T. A. Kaplan Lyons and Kaplan (1960).

II.2 Criteria for ground state degeneracies on bipartite lattices

While for Bravais lattices any ground state found via the Luttinger-Tisza method satisfies the strong constraint, this is generally not true in the case of lattices with more than one site per unit cell. In this section, we will analyze the Luttinger-Tisza solution in the special case of bipartite lattices, which consist of two sites per unit cell, i.e., they can be seen as two interpenetrating Bravais lattices with $\alpha\in\{1,2\}$. Particularly, we will show that for equivalent sublattices this solution is still exact. We will then derive a criterion for the coupling parameters such that these solutions exhibit subextensive spiral degeneracies. Our arguments here are formulated in a general way, without specifying the precise lattice geometry. In the next section, we will illustrate these considerations based on two concrete examples.

For any bipartite lattice, frustration is only possible if interactions beyond nearest-neighbor spins are present. Otherwise, neighboring spins from different sublattices can always be aligned or anti-aligned, leading to a ferromagnetic or Néel ordered state, respectively. For simplicity, in this section we restrict ourselves to first- and second-neighbor interactions only; the general argument is, however, also valid for further-neighbor couplings. Consider an $n$-dimensional lattice with a set of primitive vectors $\{\bm{a}_{i}\}$ and a basis $\bm{h}$ defining the relative shift of the two sublattices. For a bipartite lattice, the matrix $\underline{\bm{J}}(\bm{q})$ is of $2\times 2$ form. Its diagonal elements contain all the interactions between spins that are connected by primitive lattice vectors and, as a result of the inversion symmetry of Bravais lattices, are real,

[TABLE]

The off-diagonal terms ($\alpha\neq\beta$) are possibly complex,

[TABLE]

where the vectors $\bm{R}^{\text{nn}}_{\alpha\beta}$ point from an arbitrary, but fixed site on sublattice $\alpha$ to all its nearest neighbors within the other sublattice $\beta$. We now make the assumption that the two sublattices are equivalent ($J_{2}^{\alpha=1}=J_{2}^{\alpha=2}$) which means that the diagonal elements of the hermitian matrix $\underline{\bm{J}}(\bm{q})$ are equal,

[TABLE]

Diagonalization of this matrix leads to the eigenvectors and eigenvalues

[TABLE]

Here, $\phi$ is the angle between the two spins in the same unit cell given by

[TABLE]

From Eq. (13) it follows that the strong constraint can be fulfilled exactly in the case where both components of the Luttinger-Tisza eigenmode $\bm{u}_{\nu_{LT}}(\bm{q}_{LT})$ have the same absolute value. Thus, the spin configuration obtained from minimizing the eigenvalues in Eq. (17) satisfies the strong constraint for all bipartite lattices with equivalent sublattices such that Luttinger-Tisza becomes exact (note that this property is lost for inequivalent sublattices). As in the case of Bravais lattices, the solution has the form of coplanar spirals,

[TABLE]

where we use the convention $\phi_{1}\equiv\phi$ and $\phi_{2}=0$, i.e. $\phi_{1}$ is the angle between spins in the same unit cell separated by the basis $\bm{h}$.

The remaining question which we address in this subsection is under which conditions the exact solution that follows from Eq. (17) exhibits a degenerate spiral ground state manifold. For that purpose, we perform a minimization of the eigenvalues in Eq. (17) by requiring that the gradient $\nabla_{\mu}=\frac{\partial}{\partial q_{\mu}}$ ($\mu=x,y,z$) of the smallest eigenvalue must vanish,

[TABLE]

The diagonal matrix elements are of the form $J_{2}^{\star}f(\bm{q})$, where $f(\bm{q})$ is a real function containing a sum of cosine terms as in Eq. (14). The star in $J_{2}^{\star}$ (not to be confused with complex conjugation marked by an asterisk ‘∗’) indicates that this is a generalized second-neighbor coupling which may also contain longer-range couplings within the same sublattice such that the sum in Eq. (14) also includes lattice vectors between further neighbors. In the next section we will make this generalization more explicit. The absolute value of the off-diagonal elements $|J_{12}(\bm{q})|$ can be further specified using Eq. (15),

[TABLE]

For nearest neighbors, each distance vector $\bm{R}^{\text{nn}}_{12}$ in this expression consists of a sum of lattice vectors $\bm{a}_{i}$ and the basis vector $\bm{h}$, where the exponentials $e^{i\bm{q}\bm{h}}$ cancel out when being multiplied by their complex conjugate. Therefore, Eq. (21) can be recast in the more convenient form

[TABLE]

where $z$ is the coordination number of the lattice which results from contributions $e^{i\bm{q}\bm{R}^{\text{nn}}_{12}}e^{-i\bm{q}\bm{R}^{\text{nn}}_{12}}=1$. Furthermore, $pg(\bm{q})$ is a sum of cosine terms of the form $\cos\left[\bm{q}(\bm{a}_{i}-\bm{a}_{j})\right]$ where the prefactor $p$ is defined such that the smallest coefficient of cosines in $g(\bm{q})$ is one (see Sec. III for explicit examples). Using Eq. (20), possible Luttinger-Tisza eigenvalues have to fulfill the condition

[TABLE]

In the generic case, this equation corresponds to one condition for each component of $\bm{q}$ such that Eq. (23) is only fulfilled for a single spiral state. However, a special situation emerges when $f(\bm{q})=g(\bm{q})$, where the equation is either fulfilled by $\nabla f(\bm{q})=0$ or another set of solutions which is characterized by only one equation

[TABLE]

or equivalently

[TABLE]

Particularly, this single condition allows for a continuous manifold of solutions $\bm{q}$ forming closed contours (surfaces) for two (three) dimensional systems, hence, yielding spiral spin liquids. This is exactly the case for the $J_{1}$-$J_{2}$ model on the diamond and honeycomb lattices as studied in Refs. Bergman et al., 2007; Mulder et al., 2010. For the bcc lattice Attig and Trebst (2017) the condition $g(\bm{q})=f(\bm{q})$ requires the incorporation of third and fourth neighbor couplings $J_{3}$, $J_{4}$ into $J_{2}^{\star}$. All these models have in common that for $J_{1}<0$ ($J_{1}>0$) and $J_{2}^{\star}=0$ the system is first in a ferromagnetic (Néel) ground state, i.e., there is only one Luttinger-Tisza solution. However, when $J_{2}^{\star}/|J_{1}|$ is tuned beyond a certain critical ratio, Eq. (25) yields degenerate solutions. Since $f(\bm{q})$ contains only cosine terms, its leading terms are quadratic in $q_{\mu}$. Therefore, when the spiral surface with small $|\bm{q}|$ just emerges from ferromagnetic or Néel order its shape is determined by the equation of an ellipsoid $\frac{q_{x}^{2}}{c_{x}}+\frac{q_{y}^{2}}{c_{y}}+\frac{q_{z}^{2}}{c_{z}}=1$ or a sphere, as for the aforementioned models. While the above considerations hold for both signs of $J_{1}$, we will only treat the case of antiferromagnetic $J_{1}>0$ below.

It must, furthermore, be stressed that there are other ways in which Eq. (23) can result in a degenerate manifold, e.g., if two or more equations are linearly dependent. In addition, the solutions to Eq. (25) do not pose a necessary condition for absolute minima of the eigenvalues, so it has to be checked for each individual problem that this is indeed the case.

III Further models with degenerate spiral surfaces

Based on the condition presented in Sec. II.2, it is now possible to construct new models that exhibit spiral spin-liquid phases. The general recipe which works for all bipartite lattices amounts to first expressing $|J_{12}(\bm{q})|$ as in Eq. (22) and then finding the set of couplings $J_{2}^{\star}$ within the same sublattice such that $f(\bm{q})=g(\bm{q})$ is satisfied. Below we demonstrate this procedure for the square and the hcp lattice.

III.1 Square lattice

The square lattice can be decomposed into two interpenetrating square lattices which are rotated by $45^{\circ}$ and stretched by a factor $\sqrt{2}$ as compared to the original one, see Fig. 1(a). The off-diagonal elements of $\underline{\bm{J}}(\bm{q})$ take the simple form

[TABLE]

and we assume $J_{1}>0$. Because the square lattice is a Bravais lattice, $J_{12}(\bm{q})$ is real. Using the notation $q_{a}\equiv\bm{q}\bm{a}$ and $q_{b}\equiv\bm{q}\bm{b}$ where $\bm{a}$ and $\bm{b}$ are the primitive lattice vectors, the absolute value $|J_{12}(\bm{q})|$ can be written as

[TABLE]

According to Eq. (14) the diagonal elements of $\underline{\bm{J}}(\bm{q})$ have the form $J_{11}=J_{2}^{\star}f(\bm{q})$, where $J_{2}^{\star}$ contains second and (possibly) further-neighbor couplings and $f(\bm{q})$ is a sum of the cosine terms corresponding to the Fourier transforms of these couplings. The strengths of the further-neighbor couplings need to be adjusted such that the condition $f(\bm{q})=g(\bm{q})$ from Sec. II.2 is satisfied. As can be seen from the function $g(\bm{q})$ in square brackets of Eq. (27), second-neighbor couplings $J_{2}$ [generating terms $\cos(q_{a/b})$] and third-neighbor couplings $J_{3}$ [generating terms $\cos(q_{a}\pm q_{b})$] are sufficient to fulfill this condition. Furthermore, according to the prefactors of the cosine terms, the contributions from the second-neighbor couplings need to be twice as large as those from the third-neighbor couplings. In other words, $J_{2}^{\star}$ comprises a second-neighbor coupling $J_{2}$ of the strength $2J_{2}^{\star}$ and a third-neighbor coupling $J_{3}$ of the strength $J_{2}^{\star}$, as shown in Fig. 1(a).

The degenerate ground states are spirals that satisfy Eq. (25) with $z=4$ and $p=2$, i.e.,

[TABLE]

The evolution of the degenerate spiral contour as a function of the coupling ratio $g=J_{2}^{\star}/J_{1}$ is illustrated in Fig. 1(b). For $g\leq g_{\text{c}}=1/8$ there is no solution to Eq. (28), i.e. the Luttinger-Tisza wave vector is determined by the condition $\nabla f(\bm{q})=0$ and the system resides in a Néel ordered phase. Above the critical value $g_{\text{c}}=1/8$ spiral contours begin to form out of the Néel order. Close to $g_{\text{c}}$ the contours first have a circular shape but with increasing $g$ they continuously deform into the square shape of the first Brillouin zone. In the limit $g\rightarrow\infty$ the two sublattices decouple where each sublattice individually forms a square lattice model with first and second-neighbor interactions where the latter are exactly half as strong as the former. This is the well-known case of the $J_{1}$-$J_{2}$ square lattice model at the classical transition point between Néel order and collinear orderChandra and Doucot (1988). It is worth noting that the evolution of the spiral contours is very similar to the three-dimensional bcc latticeAttig and Trebst (2017), where, likewise, the degenerate manifold adopts the shape of the first Brillouin zone in the limit $g\rightarrow\infty$. One may also embed the spiral system studied here into the classical phase diagram of the $J_{1}$-$J_{2}$-$J_{3}$ model with unrestricted couplings. The line cut in parameter space defined by $g$ then corresponds to a phase boundary between a $(q,q)$ spiral and a $(q,\pi)$ spiral, both of which do not exhibit any extensive degeneracies.Moreo et al. (1990); Chubukov (1991)

III.2 hcp lattice

We finally discuss another spiral spin-liquid phase which is based on the three-dimensional hcp lattice. The hexagonal unit cell of the hcp lattice is spanned by the vectors $\bm{a}=(\frac{1}{2},-\frac{\sqrt{3}}{2},0)$, $\bm{b}=(\frac{1}{2},\frac{\sqrt{3}}{2},0)$, $\bm{c}=(0,0,c)$ and the basis $\bm{h}=(0,-\frac{\sqrt{3}}{4},-\frac{c}{2})$. One may view the hcp lattice as an ‘abab’ stacking of equilateral two-dimensional triangular lattices where the alternating sequence of these layers yields a bipartite geometry, see Fig. 2.

The off-diagonal element of the coupling matrix $\underline{\bm{J}}(\bm{q})$ is given by

[TABLE]

and its absolute value has the form

[TABLE]

with

[TABLE]

where we set $J_{1}$ to be antiferromagnetic and $q_{a}=\bm{q}\bm{a}$, $q_{b}=\bm{q}\bm{b}$, $q_{c}=\bm{q}\bm{c}$. The condition $f(\bm{q})=g(\bm{q})$ again determines the set of couplings $J_{2}^{\star}$ which is needed to generate a spiral degeneracy. By simple bookkeeping of terms one finds that one in-plane coupling of size $2J_{2}^{\star}$ as well as two couplings $J_{2}^{\star}$ and $3J_{2}^{\star}$ connecting triangular layers separated by the distance $c$ are required, see Fig. 2. The spiral surface then contains all spiral states with wave vectors $\bm{q}_{LT}$ that satisfy the equation

[TABLE]

As before, as a function of the coupling ratio $g=J_{2}^{\star}/J_{1}$ the system first shows Néel order. Above the critical coupling $g_{\text{c}}=1/12$ the condition in Eq. (32) has finite solutions such that spiral surfaces appear. The shape of the spiral surface emerging for $g>g_{\text{c}}$ also depends on the length $c=|\bm{c}|$. Interestingly, for the closest packing $c=\sqrt{8/3}$, it is not a sphere but an ellipsoid. Changing the length $c$ (which is equivalent to rescaling $q_{c}$) results in a stretching/compression of the surface in the vertical direction. By expanding $f(\bm{q})$ in Eq. (32) up to second order in $\bm{q}$ one finds that for $c=\sqrt{2/3}$ the surface for $g$ just above $g_{c}$ is exactly spherical. In Fig. 3 we plot the degenerate spiral manifolds for various values of $g$. As can be seen, for $g>0.25$ the surface cuts through the edges of the first Brillouin zone and in the limit $g\rightarrow\infty$ the surface only consists of flat planes at $\bm{q}_{z}=\pm\pi$ connected by nodal lines along the vertical edges of the first Brillouin zone.

IV Effect of quantum fluctuations

In this section, we discuss the fate of the spiral spin liquids identified in the last section when quantum fluctuations are included, i.e., when tuning the magnitude $S$ of the spins away from the classical limit. For this purpose, we employ the pseudofermion functional renormalization group (PFFRG) method for spin systems in its implementation for unrestricted $S$.Reuther and Wölfle (2010); Baez and Reuther (2017) Particularly, we calculate the magnetic correlations for the square and hcp lattice models as a function of the coupling ratio $g$ and the spin length $S$ and determine whether or not the systems develop magnetic long-range order. The results discussed below are based on the imaginary-time one-loop plus Katanin PFFRG approach for spin systems at $T=0$ as described in Ref. Baez and Reuther, 2017. The method uses an Abrikosov decomposition of each spin operator into two pseudofermions Abrikosov (1965); Abrikosov and Migdal (1970) which allows for an application of fermionic diagram techniques. The purely quartic fermionic Hamiltonians which arise after the Abrikosov decomposition of the Heisenberg interactions are notoriously difficult to solve. The PFFRG handles this problem by first introducing an artificial cutoff $\Lambda$ (the so-called RG scale) which sets the bare fermionic single-particle Green’s function to zero for Matsubara frequencies $|\omega|<\Lambda$. With this manipulation, the system’s $n$-particle vertex functions can be calculated from an infinite set of coupled integro-differential equationsWetterich (1993) which needs to be truncated for numerical solubility. Here, we use a well-established truncation scheme which takes into account one-loop and Katanin-type diagrammatic contributions.Katanin (2004) The resulting set of coupled equations is solved with the initial conditions of the vertex functions at $\Lambda\rightarrow\infty$ given by the bare interactions. Physical vertex functions are then obtained in the limit $\Lambda\rightarrow 0$ by successively integrating out energy degrees of freedom. In essence, the PFFRG performs diagrammatic resummations which systematically incorporate the large-$N$ and large-$S$ limits [$N$ determines the spins’ symmetry group $SU(N)$]. In particular, it has been proven that the leading diagrammatic contributions in $1/N$ Buessen et al. (2018b); Roscher et al. (2018) and in $1/S$ are both exactly summed up, where in the large $S$ limit the PFFRG becomes identical to the Luttinger-Tisza method.Baez and Reuther (2017) For more details on the method, please refer to Refs. Reuther and Wölfle, 2010; Baez and Reuther, 2017.

We use the PFFRG to calculate the system’s isotropic static spin correlations

[TABLE]

where $\hat{S}_{i}^{\mu}(\tau)=e^{\tau\hat{H}}\hat{S}_{i}^{\mu}e^{-\tau\hat{H}}$ and the bracket $\langle\dots\rangle_{\Lambda}$ denotes that the expectation value is computed at the RG scale $\Lambda$. Fourier-transforming $\chi_{ij}$ into momentum space, a smooth flow of the maximal $\bm{q}$-component of the susceptibility towards the physical limit $\Lambda\rightarrow 0$ indicates that the system is in a paramagnetic phase, hence, possibly realizing a quantum spin liquid. The onset of magnetic order, on the other hand, is accompanied by a divergence or a kink in the susceptibility flow. These two distinct behaviors are explained by the fact that our PFFRG formalism is invariant under a global $SU(2)$ spin-rotation. Once this symmetry is spontaneously broken due to the onset of magnetic order, the algorithm is unable to correctly describe the system’s correlations and the RG flow becomes unphysical. A true divergence at a finite critical RG scale $\Lambda_{\mathrm{c}}$ would be expected if correlations between infinitely separated lattice sites and continuous Matsubara frequencies were included. Since we use discretized frequencies and truncate the extent of the spin correlations in real space, such divergences are typically regularized to a kink in the susceptibility flow. Hence, in the case of a paramagnetic phase we can compute the ${\bm{q}}$-dependent susceptibility down to the physical limit $\Lambda\rightarrow 0$. Alternatively, from a kink or a divergence in the susceptibility flow we can infer that the system is magnetically ordered where the specific type of order is determined by the position of the magnetic Bragg peaks in $\bm{q}$ space. For the models considered here, we use a mesh of $70$ discretized points for all three Matsubara frequency arguments of the two-particle vertex.

It is worth highlighting that in the absence of an instability feature during the RG flow, our analysis below only allows us to conclude that the system is non-magnetic but does not yet imply a quantum spin-liquid phase. This is because a non-magnetic ground state may still exhibit some type of hidden order such as, e.g., dimer crystal formation. Extended PFFRG schemesReuther and Wölfle (2010); Iqbal et al. (2019, 2016b) (which will not be applied here) also allow one to probe the system with respect to dimer order. If such type of order is absent as well, a recently developed PFFRG approach is applicable which has shown some success in determining the system’s low-energy spinon band structure of the putative quantum spin liquid.Hering et al. (2019)

IV.1 Square lattice

We now present the results of our PFFRG calculations for the square lattice Heisenberg model with couplings as derived in Sec. III.1. The lattice conventions are the same as above, i.e., the primitive lattice vectors $\bm{a}$, $\bm{b}$ are defined as in Fig. 1 and have unit length $|\bm{a}|=|\bm{b}|=1$. While the PFFRG, in principle, treats an infinite lattice, finite-size effects still enter since spin correlations are only taken into account up to $10$ nearest-neighbor distances in both spatial directions and are treated as zero beyond this length. This means that, in total, each spin is coupled to $440$ surrounding spins. The free parameters of our simulations are $g=J_{2}^{\star}/J_{1}$ and $S=1/2,\,1,\,3/2,\,\dots$.

In Fig. 4 we show the momentum-resolved susceptibility $\chi(\bm{q})$ for varying $g$ and fixed $S=1/2$ as well as for fixed $g=0.25$ and $g=0.4$ and varying $S$. For a better comparison with the results from the previous sections, the Fourier transforms from $\chi_{ij}$ to $\chi(\bm{q})$ are performed with respect to one sublattice only, i.e., the Fourier sums only run over $i$, $j$ which belong to the same sublattice. As a consequence, antiferromagnetic Néel order manifests in a magnetic Bragg peak at $\bm{q}=0$. Explicitly taking into account both sublattices would result in an additional modulation of the susceptibility with a function that is periodic in the second Brillouin zone. This however, would only complicate the comparison with the classical Luttinger-Tisza results but would not yield relevant new insights.

The susceptibility plots in Fig. 4 demonstrate that even in the quantum limit $S=1/2$ the momentum distribution of the response still roughly follows the classical spiral contours. The characteristic ring-like pattern, though, sets in at slightly larger $g$ as compared to the classical limit, such that there is a small regime around $g=0.15$ where the classical system has a spiral degeneracy but the quantum system still shows a single (but broadened) peak at $\bm{q}=0$. From the flow of the susceptibility we infer that for $S=1/2$ a paramagnetic phase is realized for $g\gtrsim 0.2$ and that the system orders magnetically for smaller values of $g$ [see Fig. 5(a) which illustrates $\Lambda$-dependent susceptibilities for selected $g$ and $S$ indicating either a magnetically ordered or a non-magnetic RG flow]. Particularly, the onset of the non-magnetic phase is seen to coincide with the appearance of ring-like features in the susceptibility, revealing a close connection between both phenomena. The system remains non-magnetic in the whole parameter regime up to $g\rightarrow\infty$ where the spins on the two sublattices decouple, each realizing a $J_{1}$-$J_{2}$ square lattice model with $J_{2}=0.5J_{1}$. Previous PFFRG studies found a similar extent of the non-magnetic phase at $S=1/2$Reuther et al. (2011b) (including the limit $g\rightarrow\infty$Reuther and Wölfle (2010)), however, their focus was not on spiral properties. Furthermore, the observation that quantum fluctuations enlarge the regime of antiferromagnetic Néel correlations at the expense of spiral correlations has already been made for other systems such as Heisenberg models on the honeycomb or bcc lattices.Reuther et al. (2011a); Ghosh et al. (2019b)

A closer inspection of the susceptibility profiles for $S=1/2$ in Fig. 4 indicates some distinct differences as compared to the classical spiral contours. For $g\lesssim 0.25$ the rims in the susceptibility are found to be contracted towards the $\bm{q}=0$ point while for $g\gtrsim 0.5$ the weight is shifted more towards $\bm{q}=(\pi,\pi)$. Furthermore, the contours of strong response show a selection of dominant wave vectors due to quantum fluctuations. These effects are most pronounced around $g=0.4$ where incommensurate momenta of the form $\bm{q}=(\pm q,\pm q)$ are selected and for $g\rightarrow\infty$ where the point $\bm{q}=(\pi,\pi)$ is preferred.

Increasing the spin magnitude, there is a small range of couplings $g\simeq 0.3,\,\dots,\,0.5$ for $S=1$ where the system possibly still resides in a paramagnetic phase (the behavior of the RG flow in this regime is, though, less conclusive as compared to $S=1/2$). For all other parameter values and for higher spins, our calculations indicate a magnetically ordered ground state. As shown in the lower panels of Fig. 4, both quantum effects, i.e., contraction of spiral contours and selection of dominant wave vectors, are suppressed for increasing $S$ and the exact Luttinger-Tisza result is reproduced correctly as expected. Baez and Reuther (2017)

IV.2 hcp lattice

Next, we investigate the hcp lattice model via PFFRG. The coupling constants ranging up to fourth neighbors are defined in Sec. III.2 and our tuning parameters are again $g=J_{2}^{\star}/J_{1}$ and the spin magnitude $S$. Our calculations take into account spin correlations with a maximal length of $8$ nearest-neighbor distances such that each spin is coupled to $932$ surrounding spins. The momentum-resolved susceptibilities $\chi(\bm{q})$ are shown in Fig. 6 for varying $g$ at $S=1/2$ (top panel) and varying $S$ for $g=0.15$ and $g=0.25$ (two bottom panels). We only plot $\chi(\bm{q})$ in the plane with $q_{z}=0$ which contains most of the quantum effects we wish to discuss. Furthermore, as in Sec. IV.1, the susceptibility has been obtained by restricting to one sublattice only.

Despite the different lattice geometry and exchange couplings of this system, our observations are similar to those for the square lattice model in the previous subsection, i.e., there is an overall good agreement between the classical spiral surfaces and the regions of strong response in the calculated susceptibility. Starting with $S=1/2$, the onset of a ring-like susceptibility again occurs at a somewhat larger ratio $g$ than in the classical case. The flow of the susceptibility, which we plot in Fig. 5(b) for selected parameter values, implies that there is an extended range $g\simeq 0.2,\,\dots,\,1$ where the spin-$1/2$ system is in a paramagnetic phase. Interestingly, in contrast to the square lattice model of the previous subsection, here we also identify a small parameter regime around $g=0.15$ where spiral surfaces are clearly visible in our PFFRG results but the system is still magnetically ordered. Within the paramagnetic phase the spiral contours in the $q_{z}=0$ plane first show a pattern of kagome-like streaks at $g=0.25$ which then contract towards the corners of the first Brillouin zone. We again observe small deviations between the classical degeneracies and the susceptibility distribution in the quantum case. Particularly, for small $g$ (see, e.g., $g=0.15$) the rings in the susceptibility are smaller than classically expected while for larger $g\gtrsim 0.375$ they quickly become pinned to the corners of the first Brillouin zone. Selection effects at $S=1/2$ seem less pronounced compared to the square lattice model, nevertheless, a certain selection takes place outside the $q_{z}=0$ plane which is not visible in our plots. In the limit $g\rightarrow\infty$, the hcp layers decouple into an ‘aaa’ stacking of triangular layers with weak indication of magnetic order.

With increasing $S$ the non-magnetic phase quickly shrinks, e.g., at $S=1$, there only exists a small regime $g\simeq 0.2,\,\dots,\,0.4$ where the system is possibly non-magnetic. For $S>1$ only magnetic phases remain and the susceptibility distribution shows increasingly better agreement with the Luttinger-Tisza result. However, we also emphasize that three-dimensional systems with large spin magnitudes as studied here are rather prone to finite-size effects within the PFFRG. This is because spin correlations have a strong tendency to become long-range when quantum fluctuations die out for $S\rightarrow\infty$. Furthermore, for three-dimensional systems we have to reduce the distance of the longest spin correlations in our algorithm to become numerically feasible. As a consequence, there may easily arise a situation where the correlation length is much larger than the distance of the longest correlations taken into account. This effect is most pronounced in our susceptibility plot for $g=0.25$ and $S=5$ showing a spurious selection of wave vectors at the midpoint positions of the edges of the Brillouin zone which are not expected at large $S$.

V Conclusions

In summary, this work investigates the capability of classical spiral spin liquids to induce quantum spin-liquid behavior when the spin magnitude is decreased down to $S=1/2$. For several spiral spin-liquid systems such as classical Heisenberg models on the honeycomb, diamond, and bcc lattices the onset of a non-magnetic quantum phase at $S=1/2$ has already been proposed in various numerical studies.Buessen et al. (2018a); Reuther et al. (2011a); Zhu et al. (2013); Gong et al. (2013); Zhang and Lamas (2013); Li et al. (2012); Ghosh et al. (2019b) However, due to the lack of a common principle behind the precise arrangements of couplings in these systems, the family of classical spin models yielding spiral degeneracies has been very small so far. To overcome this deficiency, we have first developed a general and exact procedure to construct new spiral spin-liquid phases on bipartite lattices. As an example, we have demonstrated this approach based on the two-dimensional square lattice and the three-dimensional hcp lattice, both of which exhibit a one-parameter family of exchange couplings where continuous spiral degeneracies exist. We have further studied the impact of quantum fluctuations on these systems by employing the PFFRG method. We find that large portions of the classical spiral spin-liquid regime become non-magnetic quantum phases when the spin length is reduced to $S=1/2$. However, even in this extreme quantum limit spiral correlations still determine the system’s magnetic properties on short length scales, as evidenced by pronounced rims of strong response at approximately the same $\bm{q}$-space locations as the classical spiral contours. Depending on the precise coupling ratios, the susceptibility is not evenly distributed along these rims but shows small maxima which indicate a quantum order-by-disorder effect.

While the spin models investigated in this work are interesting candidate systems for stabilizing exotic non-magnetic quantum phases, it needs to be emphasized again that based on our current analysis we cannot draw any definite conclusion about whether or not they indeed realize a quantum spin liquid. A well-known alternative is the occurrence of spontaneous dimer ordering which may also be checked within a modified PFFRG approach.Reuther and Wölfle (2010); Iqbal et al. (2019, 2016b) Such an analysis, however, is beyond the scope of the current work. Eventually, it would be desirable to realize the above lattice geometries and sets of exchange couplings in a real material. As a complicating fact, though, both models exhibit longer-range exchange couplings beyond second neighbors which need to appear in fixed ratios to each other. While it might be difficult to exactly realize these ratios in a material, our procedure for designing spiral spin liquids is very general and may be applied to all bipartite lattices. Particularly, as a future direction of research, it would be interesting to search for new spiral spin-liquid models more systematically and identify simpler and more realistic arrangements of exchange couplings.

Acknowledgements

We would like to thank Stef Koenis for his feedback and helpful discussions. The PFFRG simulations were performed on the tron cluster at Freie Universität Berlin. This work was partially supported by the German Research Foundation within the CRC 183 (project A02).