Phase diagram of a distorted kagome antiferromagnet and application to Y-kapellasite

Abstract

We investigate the magnetism of a previously unexplored distorted spin-1/2 kagome model consisting of three symmetry-inequivalent nearest-neighbor antiferromagnetic Heisenberg couplings J_⬡, J, and $J^{\prime}$, and uncover a rich ground state phase diagram even at the classical level. Using analytical arguments and numerical techniques we identify a collinear $\overrightarrow{Q}=0$ magnetic phase, two unusual non-collinear coplanar $\overrightarrow{Q}=(1/3,1/3)$ phases and a classical spin liquid phase with a degenerate manifold of non-coplanar ground states, resembling the jammed spin liquid phase found in the context of a bond-disordered kagome antiferromagnet. We further show with density functional theory calculations that the recently synthesized Y-kapellasite Y₃Cu₉(OH)₁₉Cl₈ is a realization of this model and predict its ground state to lie in the region of $\overrightarrow{Q}=(1/3,1/3)$ order, which remains stable even after the inclusion of quantum fluctuation effects within variational Monte Carlo and pseudofermion functional renormalization group. The presented model opens a new direction in the study of kagome antiferromagnets.

Introduction

The kagome lattice is arguably one of the most important two-dimensional (2D) lattices for the study of magnetic frustration. It is characterized by a complex phase diagram including magnetically ordered regimes and proposed quantum spin liquid phases¹, has rich magnetization dynamics², and supports some of the best-studied quantum spin liquid candidates like herbertsmithite ZnCu₃(OH)₆Cl₂^3,4,5,6. In a more technical context, the study of the antiferromagnetic Heisenberg model on the kagome lattice has been a fertile ground for the development and benchmarking of theoretical methods. Notably, the competition between density matrix renormalization group^7,8,9, variational Monte Carlo (VMC)^10,11, and tensor networks¹² type methods, with the aim of resolving the nature of the spin liquids supported by the kagome lattice, has been a fervent area of research for many years.

All these intense research activities have mainly focused on the ideal kagome structure. In contrast, distortions of this lattice have been studied much less, even though they are realized in some magnetic compounds, and their physical phenomenology may be even richer than for the standard kagome lattice. In some cases, like volborthite Cu₃V₂O₇(OH)₂·2H₂O^13,14, the distortion leads to a new 2D lattice which is still highly frustrated and possibly has a spin liquid ground state¹⁵. In Rb₂Cu₃SnF₁₂, the deformed kagome lattice leads to a pinwheel valence bond solid¹⁶. Other kinds of distortions lower the rotational symmetry of the lattice and lead to kagome strips^17,18. Even the low-temperature structure of herbertsmithite bears some signatures of distortion^19,20,21.

The focus of the present work lies on an unusual and previously unexplored distortion of the kagome lattice which is realized in the recently synthesized variant of herbertsmithite, namely Y₃Cu₉(OH)₁₉Cl₈. The distorted lattice structure consists of three symmetry-inequivalent nearest-neighbor kagome bonds forming a nine-site unit cell. Analyzing the corresponding Heisenberg model as a function of its two coupling ratios, using analytical arguments and numerical techniques, we find a surprisingly rich ground state phase diagram, even at the classical level. A first notable observation is that large parts of the phase diagram represent an unusual coplanar spin state with a commensurate magnetic wave vector $\overrightarrow{Q}=(1/3,1/3)$. This type of ordered state requires a 27 atom magnetic unit cell. Furthermore, in an extended regime around the standard undistorted kagome lattice, an even more complex classical spin liquid phase is identified, which cannot be characterized by any specific wave vector. It bears similarities with the well-known classical spin liquid on the undistorted kagome lattice in the sense that its low energy states follow from a set of spin constraints for each triangle²². In contrast, however, the ground state in the distorted case is found to be generally non-coplanar and, hence, resembles the jammed spin liquid investigated in ref. ²³.

After discussing in detail the general magnetic phenomena of this distorted kagome lattice, the second focus of this paper is on the specific case of Y₃Cu₉(OH)₁₉Cl₈ where this model is likely realized. This material was discovered in an attempt of electron doping the Cu 3d states in herbertsmithite with the aim of placing the Fermi level at the symmetry-protected Dirac crossing of the Cu d-bands in the kagome lattice^24,25,26. In herbertsmithite-type copper hydroxy-halides, however, the larger charge provided by Y³⁺ compared to Zn²⁺ is always compensated by the incorporation of additional hydroxide or halide ions, preserving the antiferromagnetic insulator nature of the Cu²⁺ layers. Even if charge doping remains elusive in these systems, the newly discovered by-product in form of Y₃Cu₉(OH)₁₉Cl₈ appears to be of great interest in and of itself.

In Y₃Cu₉(OH)₁₉Cl₈, the Y³⁺ ions are placed in the center of the hexagon of the kagome lattice, making it a material that is structurally similar to kapellasite^27,28,29, haydeeite³⁰, or centennialite^31,32. We, therefore, name the system Y-kapellasite. Note that there is a closely related compound YCu₃(OH)₆Cl₃ with an ideal kagome lattice but disordered Y positions³³. The latter orders at T_N = 15 K in a $\overrightarrow{Q}=0$ structure with negative spin chirality^34,35 which has been attributed to a strong Dzyaloshinskii–Moriya (DM) interaction³⁶. In contrast, Y-kapellasite remains dynamical down to much lower temperatures than YCu₃(OH)₆Cl₃³⁷; it has a broad feature at T = 2 K in the specific heat but muon spin resonance (μSR) on powder samples seems to indicate the absence of any static magnetic order, although disorder effects may play a role. Recently, the coexistence of magnetic order and persistent spin dynamics has been suggested for different samples of Y-kapellasite³⁸.

By extracting the Heisenberg Hamiltonian of Y-kapellasite using total energy mapping from density functional theory (DFT) calculations, we find that the three couplings on the symmetry-inequivalent nearest-neighbor kagome bonds dominate, with negligible longer range interactions. We may, hence, place Y-kapellasite in the region of $\overrightarrow{Q}=(1/3,1/3)$ order in the classical ground state phase diagram obtained here. Investigating the corresponding spin-1/2 model within VMC and pseudofermion functional renormalization group (PFFRG) we argue that quantum fluctuations are not sufficiently strong to suppress the long-range magnetic order. Accordingly, our semiclassical spin-wave analysis provides a realistic approximation of the system’s excitation spectrum which will be useful for comparison with future experimental data.

Results

Spin Hamiltonian

The model investigated in this work is a variant of the standard nearest-neighbor kagome Heisenberg model, but with three distinct nearest-neighbor couplings, which we call J, J_⬡, and $J^{\prime}$ [see Fig. 1a]. We will later argue that this model approximates well the microscopic interactions in Y-kapellasite. The Heisenberg Hamiltonian can be written as

$${{{\mathcal{H}}}}=\mathop{\sum}\limits_{\langle i,j\rangle }{J}_{ij}{\overrightarrow{S}}_{i}\cdot {\overrightarrow{S}}_{j}\ ,$$

(1)

where ${\overrightarrow{S}}_{i}$ are the spin degrees of freedom (which, below, are either chosen as spin-1/2 operators or as classical normalized vectors) and J_ij is given by J, J_⬡, or $J^{\prime}$, depending on the bond. All these couplings are assumed to be positive (antiferromagnetic). It is clear that $J={J}_{\hexagon}=J^{\prime}$ leads back to the standard undistorted nearest-neighbor kagome model. As a consequence of the broken translational symmetry of the kagome lattice, the system’s periodic structure is described by a decorated triangular lattice with a unit cell of nine sites [see Fig. 1a]. We can distinguish two inequivalent sets of sites inside the unit cell, which are not connected by point group symmetries and form two distinct sublattices: sublattice A is made of the six sites connected by J_⬡ [the vertices of the red hexagons of Fig. 1a]; sublattice B is made of the remaining three sites. Also, note that the model is invariant under exchanging J and $J^{\prime}$ followed by a reflection with respect to the ${\overrightarrow{a}}_{1}$ axis.

Fig. 1: Distorted kagome lattice and classical $\overrightarrow{Q}=(1/3,1/3)$ magnetic order.

a Schematic illustration of the three exchange couplings characterizing the effective Heisenberg Hamiltonian for Y-kapellasite (J_⬡, J, and $J^{\prime}$) shown in red, blue, and green, respectively. The presence of three different couplings breaks the translational symmetry of the kagome lattice and leads to a decorated triangular lattice with an enlarged unit cell of nine sites, here represented by the black hexagon (Wigner–Seitz cell). The Hamiltonian of the system is periodic under translations along the Bravais vectors ${\overrightarrow{a}}_{1}$ and ${\overrightarrow{a}}_{2}$ and the sites within the unit cell are divided into sublattices A (hollow symbols) and B (solid symbols). Due to the different values of the three exchange terms, the D₆ point group symmetry of the kagome lattice is broken down to C₆. b Pictorial view of the reciprocal space. The blue arrows represent the unit vectors of the reciprocal space (${\overrightarrow{b}}_{1}$ and ${\overrightarrow{b}}_{2}$). The gray hexagons tiling the reciprocal space depict the first Brillouin zone of the lattice, while the black dashed hexagon delimits the so-called extended Brillouin zone. Some of the high symmetry points are marked with black dots. Finally, red lines represent the path along which the magnon dispersion is plotted in Fig. 3. c Classical $\overrightarrow{Q}=(1/3,1/3)$ magnetic order for $J \,>\, J^{\prime}$ (red region of Fig. 2). The orientations of the spins are fully specified by the angle ϕ between neighboring spins in the J_⬡-hexagons, as outlined in the main text (here we take the value of ϕ for the case $J^{\prime} =0$ and J_⬡ = J). In this figure, the spins are arranged in the xy-plane and their orientation is represented by the angle with respect to the S_x axis. The red, blue, and green colors of the spins of sublattice A help visualizing the $\overrightarrow{Q}=(1/3,1/3)$ pattern. d The classical $\overrightarrow{Q}=(1/3,1/3)$ magnetic order of Fig. 1c in the J ≫ J_⬡ limit ($J^{\prime} =0$). The spins form antiferromagnetic trimers along the J-bonds (depicted in black). The trimers are arranged in an effective kagome lattice structure and their orientations, highlighted by the three different colors, follow the $\sqrt{3}\times \sqrt{3}$ pattern⁶⁹.

Classical phase diagram

In Fig. 2, we summarize the classical ground-state phase diagram of the Heisenberg model on the distorted kagome lattice [Eq. (1)] as a function of the ratios J/J_⬡ and $J^{\prime} /{J}_{\hexagon}$, which has been obtained by combining analytical arguments, iterative minimization and classical Monte Carlo calculations (see Methods). At the classical level, we observe (i) a collinear $\overrightarrow{Q}=0$ magnetic phase, (ii) two non-collinear coplanar magnetic phases, both labeled as $\overrightarrow{Q}=(1/3,1/3)$ order, separated by (iii) a classical spin-liquid phase with a degenerate manifold of non-coplanar ground states, which in the context of a bond-disordered kagome antiferromagnet was dubbed a “jammed spin liquid” phase²³.

Fig. 2: Classical phase diagram of the distorted kagome model.

We note that the phase diagram is symmetric under the exchange of the axes (i.e., $J\leftrightarrow J^{\prime}$), as a consequence of the symmetry of the Hamiltonian. The $\overrightarrow{Q}=(1/3,1/3)$ magnetic order of the red region is depicted in Fig. 1c, and related to the $\overrightarrow{Q}=(1/3,1/3)$ order of the blue region by a mirror reflection with respect to ${\overrightarrow{a}}_{1}$. Inside the gray area, the system features a classical spin-liquid phase with degenerate non-coplanar ground states, as discussed in the main text. The $\overrightarrow{Q}=0$ magnetic order can be viewed as parallel spins on the same sublattice, while on different sublattices the spins are anti-aligned. We note that the axes change between the four quadrants of this plot. The empty and filled squares indicate two possible sets of couplings for Y-kapellasite.

Even without any prior knowledge of the precise nature of the different phases, there exists a simple argument that determines the location of the phase boundaries. To this end, we employ the analytical procedure of ref. ²³ where a bond disordered Heisenberg model on the kagome lattice was studied. In the first step, we rewrite the Hamiltonian of Eq. (1) as

$${{{\mathcal{H}}}}=\frac{1}{2}\mathop{\sum}\limits_{\bigtriangleup }{\left({\overrightarrow{L}}_{\bigtriangleup }\right)}^{2}+\,{{\mbox{const.}}}$$

(2)

where the sum runs over all triangles formed by nearest-neighbor bonds of the kagome lattice (both up and down triangles are considered). We define

$${\overrightarrow{L}}_{\bigtriangleup }=\sqrt{\frac{{J}_{ij}{J}_{ik}}{{J}_{jk}}}{\overrightarrow{S}}_{i}+\sqrt{\frac{{J}_{ji}{J}_{jk}}{{J}_{ik}}}{\overrightarrow{S}}_{j}+\sqrt{\frac{{J}_{ki}{J}_{kj}}{{J}_{ij}}}{\overrightarrow{S}}_{k}\ ,$$

(3)

where i, j, k ∈ Δ are the three sites forming a triangle. In our distorted model, all triangles of the lattice are formed by one J_⬡, one J, and one $J^{\prime}$ coupling [see Fig. 1a]. Thus, by an appropriate choice of the i, j, k labels of Eq. (3), we can write

$${\overrightarrow{L}}_{\bigtriangleup }=\sqrt{JJ^{\prime} /{J}_{\hexagon}}{\overrightarrow{S}}_{i}+\sqrt{J{J}_{\hexagon}/J^{\prime} }{\overrightarrow{S}}_{j}+\sqrt{J^{\prime} {J}_{\hexagon}/J}{\overrightarrow{S}}_{k}$$

(4)

for all triangles. From Eq. (2) it immediately follows that any spin configuration that fulfills the condition ${\overrightarrow{L}}_{\bigtriangleup }=0\,\,\forall \bigtriangleup$ is a ground state of the system. However, depending on the values of the couplings J, J_⬡, and $J^{\prime}$, it may occur that ${\overrightarrow{L}}_{\bigtriangleup }=0$ is impossible for any triangle when one term on the right-hand side of Eq. (4) dominates so strongly that it cannot be compensated by the other two terms.

Restricting to an isolated triangle, it is easy to show that ${\overrightarrow{L}}_{\bigtriangleup }=0$ can only be fulfilled if

$$J/{J}_{\hexagon}\ \le \ J^{\prime} /({J}_{\hexagon}-J^{\prime} )\ ,\quad J^{\prime} \ \le \ \min (J,{J}_{\hexagon})\ ,$$

(5a)

$$J/{J}_{\hexagon}\ \ge \ J^{\prime} /({J}_{\hexagon}+J^{\prime} )\ ,\quad J\ \le \ \min (J^{\prime} ,{J}_{\hexagon})\ ,$$

(5b)

$$J/{J}_{\hexagon}\ \ge \ J^{\prime} /(J^{\prime} -{J}_{\hexagon})\ ,\quad {J}_{\hexagon}\ \le \ \min (J,J^{\prime} )\ .$$

(5c)

These conditions define the phase boundaries in Fig. 2. In the regions where an isolated triangle cannot satisfy ${\overrightarrow{L}}_{\bigtriangleup }=0$, the system realizes one of the aforementioned coplanar phases [$\overrightarrow{Q}=(1/3,1/3)$ order and $\overrightarrow{Q}=0$ order]. On the other hand, in the region where an isolated triangle can fulfill Eq. (5), we observe a classical spin-liquid phase. We note that analogous phase boundaries characterize the classical phase diagram of the square-kagome antiferromagnet³⁹.

Coplanar orders

We start our discussion of the classical ground states with the coplanar phases where ${\overrightarrow{L}}_{\bigtriangleup }=0$ is necessarily violated. The rewritten Hamiltonian in Eq. (2) still implies that these phases form in a way that minimizes ${({\overrightarrow{L}}_{\bigtriangleup })}^{2}$. To simplify the investigation, we first restrict ourselves to the case $J^{\prime} =0$ where the $\overrightarrow{Q}=(1/3,1/3)$ phase is realized. In the phase diagram of Fig. 2, this corresponds to the leftmost vertical axis and it will turn out to provide a good approximation of the exchange couplings of Y-kapellasite determined by the ab initio DFT calculations (marked with squares in the figure).

In the limit $J^{\prime} =0$, the model consists of a lattice of hexagons, made of sublattice A sites, which are connected to each other through the J-trimers involving sublattice B sites [Fig. 1a]. Note that the middle spin of each trimer is fixed in the direction opposite to the sum of the edge spins. The magnetic order realized along the $J^{\prime} =0$ line is depicted schematically in Fig. 1c. The spins are coplanar and form a periodic configuration with momentum $\overrightarrow{Q}=(1/3,1/3)$ (in units of the reciprocal lattice vectors ${\overrightarrow{b}}_{1}$ and ${\overrightarrow{b}}_{2}$). This momentum corresponds to the K point of the Brillouin zone of the lattice (and the symmetry-related points), see Fig. 1b. Within a given unit cell, the spins of sublattice A form an alternating pattern around the J_⬡-hexagons: the spins on even and odd sites are ferromagnetically aligned along two different directions, which are rotated with respect to each other by an angle ϕ. The orientations of the spins on the remaining sites (i.e., sublattice B), which are only two-coordinated, are uniquely determined by the value of the angle ϕ. Thus, in the limit $J^{\prime} =0$, we can express the classical energy per site of the $\overrightarrow{Q}=(1/3,1/3)$ order as a simple function of ϕ

$$E/N=\frac{2}{3}\left[{J}_{\hexagon}\cos (\phi )+J\cos \left(\frac{\phi }{2}+\frac{\pi }{3}\right)\right]$$

(6)

Minimization yields optimal angles ϕ that go from ϕ = π in the strong hexagon limit J ≪ J_⬡ to $\phi =\frac{4\pi }{3}$ in the trimer limit J ≫ J_⬡ (see Fig. 1d and Supplementary Note 2).

Going away from the $J^{\prime} =0$ limit, the $\overrightarrow{Q}=(1/3,1/3)$ order extends for a finite region along the $J^{\prime} /{J}_{\hexagon}$ axis (red area in Fig. 2), which is bounded by the onset of the classical spin-liquid phase. Within this region, the numerical minimization of the classical energy shows that the spin pattern is unchanged with respect to the one shown in Fig. 1c and the orientation of the spins is still determined only by the angle ϕ between the spins of sublattice A. As already mentioned, the phase diagram is invariant under the exchange of J and $J^{\prime}$, and a corresponding $\overrightarrow{Q}=(1/3,1/3)$ order is observed also in the proximity of the J = 0 limit (blue area in Fig. 2). The two $\overrightarrow{Q}=(1/3,1/3)$ ordered phases, one for $J\, >\, J^{\prime}$ and the other for $J\, < \,J^{\prime}$, are transformed into each other by a mirror reflection with respect to ${\overrightarrow{a}}_{1}$. In the numerical calculations, the two phases can be distinguished by their spin susceptibility in momentum space, defined as

$${\chi }_{\vec{k}}=\frac{1}{N}\mathop{\sum}\limits_{i,j}{e}^{i\vec{k}\cdot ({\vec{r}}_{i}-{\vec{r}}_{j})}\left\langle {\overrightarrow{S}}_{i}\cdot {\overrightarrow{S}}_{j}\right\rangle \ ,$$

(7)

where ${\vec{r}}_{i}$ is the position of site i in the kagome lattice, N is the total number of sites and the brackets 〈…〉 denote an appropriate average. In the $\overrightarrow{Q}=(1/3,1/3)$ phase at $J\, > \,J^{\prime}$, ${\chi }_{\vec{k}}$ displays high intensity peaks at the ${K}_{2}^{\prime}$ points, while in the ordered phase at $J \,<\, J^{\prime}$, its maxima are located at the ${K}_{3}^{\prime}$ points [see Fig. 1b]. In the limit where ${J}_{\hexagon}\ll J,J^{\prime}$ the system is no longer frustrated since each coupled neighbor of an A site is a B site and vice versa. The ground state order is, hence, given by a simple collinear $\overrightarrow{Q}=0$ state where the two sublattices have opposite spin orientations. This regime is marked green in Fig. 2.

As it represents a previously unexplored magnetic state, it is interesting to study the classical spin wave dispersion of the $\overrightarrow{Q}=(1/3,1/3)$ order. In Fig. 3, we show the spin wave spectra for the $\overrightarrow{Q}=(1/3,1/3)$ magnetic order in three paradigmatic regimes (with $J^{\prime} =0$ in each case): J ≪ J_⬡ (strong hexagon limit) [Fig. 3a], J = J_⬡ [Fig. 3b], and J ≫ J_⬡ (strong trimer limit) [Fig. 3c]. In all three cases, the spin-wave spectrum has gapless modes at Γ and K points. For J ≪ J_⬡, J = J_⬡, we observe a finite gap between the low-lying magnon band (which gives rise to three branches when folded) and the higher bands. In the strong hexagon limit, where the system is made of weakly coupled hexagons forming a triangular pattern, the excitation spectrum at low energies resembles that of the triangular lattice antiferromagnet (see Supplementary Note 2 and Supplementary Fig. 3). The gap between the low-lying branches and the higher ones closes when the ratio J/J_⬡ is sufficiently large, as shown by the J ≫ J_⬡ case [Fig. 3c]. In this limit (strong trimer limit), the system is described by trimers of spins forming an effective kagome lattice, and the spin wave spectrum resembles that of the kagome $\sqrt{3}\times \sqrt{3}$ magnetic order (see Supplementary Note 2 and Supplementary Fig. 3).

Fig. 3: Spin wave theory results.

Calculated spin-wave dispersion within linear spin-wave theory along the (ξ, ξ)-direction [from Γ to Γ″, in Fig. 1b] for three different cases: a J/J_⬡ = 1/5, b J/J_⬡ = 1, cJ/J_⬡ = 25, where we set ${J}^{\prime}=0$ in each case. For all the choices of couplings considered here, the system is in the $\overrightarrow{Q}=(1/3,1/3)$ ordered phase. The energy scale is set by J_⬡ in panel (a), and by J in panels (b) and (c). We note that in the J ≪ J_⬡ (a) and J = J_⬡ (b) cases the lowest band (three-times folded) is separated from the higher bands by an energy gap. This gap closes upon increasing the ratio J/J_⬡. d Calculated spin-wave dispersion and intensity within linear spin wave theory along the path Γ–K–${K}_{1}^{\prime}$–Γ^″–${K}_{2}^{\prime}$–${K}_{3}^{\prime}$–Γ″ [see Fig. 1b] for the ab initio estimated Heisenberg couplings J = 154.4 K, J_⬡ = 134.2 K and ${J}^{\prime}=8.7$ K (S_XRD structure). The spectral intensity is given by the perpendicular component of the spin dynamical structure factor, which is related to the cross section of unpolarized neutron scattering experiments⁷⁰ (the overall color scale is in arbitrary units). The low-energy spectral intensity at the ${K}_{1}^{\prime}$ (${K}_{3}^{\prime}$) point is approximately 40% (20%) of the maximal intensity, located at the ${K}_{2}^{\prime}$ point. We apply a Gaussian broadening with a standard deviation of 0.2 meV.

Classical spin liquid phase

Finally, we consider the regime where the three couplings generally enable the fulfillment of ${\overrightarrow{L}}_{\bigtriangleup }=0$ in Eq. (4) and we unveil the ground state nature of this intriguing phase. Even though satisfying ${\overrightarrow{L}}_{\bigtriangleup }=0$ in an isolated triangle is possible, this does not immediately imply that achieving ${\overrightarrow{L}}_{\bigtriangleup }=0$ in each individual triangle of the full system is a trivial task. In ref. ²³ a generic bond disordered kagome system was investigated for which it was shown that ${\overrightarrow{L}}_{\bigtriangleup }=0$ can be satisfied in each triangle. Furthermore, the authors constructed global ground states where each triangle may realize up to two possible spin configurations that locally obey ${\overrightarrow{L}}_{\bigtriangleup }=0$ resulting in an extensive, but discretely degenerate classical spin liquid forming a set of the ground states with the cardinality , which was named “jammed spin liquid”.

We performed iterative minimization of the classical energy to confirm the presence of a degenerate manifold of non-coplanar ground states with ${\overrightarrow{L}}_{\bigtriangleup }=0$ within the gray region of Fig. 2 (see Methods). The main results of the minimization are summarized in Fig. 4 a, where we plot the value of $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle$ (averaged over triangles) as a function of the exchange couplings, for a finite-size cluster with N = 27 sites (see Supplementary Fig. 1, and Supplementary Fig. 3 for analogous results on a N = 36 site cluster). The optimal solutions are actually found to be identical for all triangles. Its square value yields the residual energy per triangle with respect to the ideal ground state with ${\overrightarrow{L}}_{\bigtriangleup }=0$. We observe that in most of the regions delimited by the boundaries of Eq. (5), we obtain a degenerate set of ground states with $| {\overrightarrow{L}}_{\bigtriangleup }{| }^{2}=0$ for each triangle (within machine precision). Indeed, starting the iterative minimization with different initial points, we find several independent minima which cannot be connected to each other by lattice symmetries and global spin rotations, thus confirming the large degeneracy of the classical ground state, as predicted in ref. ²³. We also note that the ${\overrightarrow{L}}_{\bigtriangleup }=0$ solutions found by the minimization are, in general, non-coplanar, and can be exploited to construct ${\overrightarrow{L}}_{\bigtriangleup }=0$ ground states for larger systems, by simply using the N = 27 sites cluster as effective unit cells for type-II clusters (see Supplementary Note 1). However, close to the boundaries of the non-coplanar region, we find a number of points in the phase diagram where achieving ${\overrightarrow{L}}_{\bigtriangleup }=0$ by numerical minimization was not possible. We are not able to provide a final statement whether the $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle\, >\, 0$ points close to the boundaries are an artifact of the finite-size calculations (and boundary conditions²³), or whether they belong to the neighboring $\overrightarrow{Q}=(1/3,1/3)$ ordered region, which may extend slightly beyond the ideal boundaries of Eq. (5). Indeed, it is worth noting that the best $\overrightarrow{Q}=(1/3,1/3)$ solution at the analytical phase boundary with the spin-liquid phase [i.e., when the equality of Eq. (5a) or (5b) holds] always has finite residual energy, i.e., $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle\, >\, 0$ (except for the extreme cases where one coupling is zero). Thus, a continuous deformation of the $\overrightarrow{Q}=(1/3,1/3)$ order to the $| {\overrightarrow{L}}_{\bigtriangleup }| =0$ spin-liquid phase cannot take place at the precise location of the analytical boundaries. Therefore, either the transition is not continuous, or the position of the phase boundaries is slightly shifted with respect to the analytical conditions of Eq. (5a) and (5b). Nevertheless, except for the precise location of the transition, our numerical minimization confirms the presence of an extended classical spin-liquid phase, characterized by degenerate non-coplanar ground states.

Fig. 4: Local ground state constraint within the classical spin liquid phase.

a Value of $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle$ (averaged over triangles) in the optimal classical ground state found by iterative minimization. The results are obtained in the classical spin liquid phase on an N = 27 sites cluster (see Supplementary Fig. 1). Most of the points fulfill the ${\overrightarrow{L}}_{\bigtriangleup }=0$ constraint for the ground state. The small number of points with $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle\, > \,0$ are found to have a coplanar non-degenerate ground state which corresponds to the nearby $\overrightarrow{Q}=(1/3,1/3)$ order. b Phase diagram from classical Monte Carlo simulations. We depict the average of $| {\overrightarrow{L}}_{\bigtriangleup }|$ for 10 simulated systems with N = 675 sites (type-II cluster, L = 5, cf. Supplementary Note 1) at the temperature $T=0.001\ {J}_{\max }$, where ${J}_{\max }=\max ({J}_{\hexagon},J,J^{\prime} )$. The $\overrightarrow{Q}=(1/3,1/3)$ ($\overrightarrow{Q}=0$) ordered phase, with maximal spin susceptibility at the ${K}_{2}^{\prime}$ or ${K}_{3}^{\prime}$ (Γ) points in the Brillouin zone, lies outside the region defined by Eqs. (5) where ${\overrightarrow{L}}_{\bigtriangleup }$ can never be zero. Inside the region where ${\overrightarrow{L}}_{\bigtriangleup }$ can potentially vanish, we still find finite but small values. The maximum value in the plot, ${\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle }_{\max }\simeq 2.83$, is found at the points of maximal distortion, i.e., $J^{\prime} \ (J)=10\ {J}_{\hexagon}=100\ J\ (J^{\prime} )$. The logarithmic color function scales as ${{\mathrm{ln}}}\,(100\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle +1)/{{\mathrm{ln}}}\,(100{\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle }_{\max }+1)$. The empty and filled squares indicate two possible sets of couplings for Y-kapellasite.

In addition to energy minimization, we performed classical Monte Carlo simulations in the low-temperature limit (see Methods), computing the value of $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle$ in the full phase diagram, as shown in Fig. 4b. The Monte Carlo results confirm the presence of a region of non-coplanar ground states within the boundaries of Eq. (5). In this region, the value of $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle$ is found to be clearly smaller than in the rest of the phase diagram, where the $\overrightarrow{Q}=(1/3,1/3)$ and $\overrightarrow{Q}=0$ orders are observed. The finite value of $\langle | {\overrightarrow{L}}_{\bigtriangleup }| \rangle$ within the non-coplanar region can be ascribed to the effect of finite temperature. To further characterize the properties of this phase, we compute the spin susceptibility [Eq. (7)]. As shown in Fig. 5f, the spin–spin correlations in the non-coplanar phase cannot be described by any particular wave vector, but rather by a distribution of wave vectors.

Fig. 5: Classical Monte Carlo results.

a–e Spin susceptibility in momentum space at different temperatures from classical Monte Carlo simulations for the S_XRD structure on a cluster with N = 7803 spins (type-II cluster with L = 17, see Supplementary Note 1), with γ = 0.001 (see Methods for details). We only consider finite real-space correlations within a circle with a radius of 50 nearest-neighbor distances around each spin. The extended Brillouin zone is depicted as a black hexagon, cf. Fig. 1b. f Spin susceptibility for J/J_⬡ = 0.5 and $J^{\prime} /{J}_{\hexagon}=0.45$, within the non-coplanar phase. The susceptibility has been computed by classical Monte Carlo calculations for ten N = 7803 sites clusters (type-II, L = 17, see Supplementary Note 1) at T = 0.001 J_⬡.

Magnetic Hamiltonian of Y-kapellasite

We now concentrate on the specific case of Y-kapellasite and investigate the magnetic properties of its spin Hamiltonian in more detail, also including the effects of quantum fluctuations. We start by performing ab initio DFT calculations to confirm that Y-kapellasite supports the spin model of Fig. 1a and determine the precise values of the coupling constants. We used both published crystal structures, the one determined by X-ray diffraction of single crystals²⁶, and the structure of Y₃Cu₉(OD)₁₉Cl₈ determined by neutron diffraction on powder samples³⁷. We consider the former structure more reliable than the latter (we follow the privately communicated assessment of P. Puphal that the single crystal samples of ref. ²⁶ are purer and less strained than the deuterated powders of ref. ³⁷), and we will refer to the single crystal structure²⁶ as S_XRD structure and to the powder structure³⁷ as S_ND structure (see Methods for more details). Note that our analysis below is valid for both structures.

The three largest couplings J, J_⬡, and ${J}^{\prime}$ (all antiferromagnetic) are shown in Fig. 6b, as obtained by energy mapping for the S_XRD structure of Y-kapellasite. The couplings are tabulated in Supplementary Table 1. These three couplings form the distorted kagome lattice illustrated in the inset of Fig. 6b. Our reasoning that the relevant exchange couplings for Y-kapellasite are just the three nearest neighbors of the distorted kagome lattice is based on extensive energy mapping for seven and, for additional confidence, 24 neighbors up to Cr-Cr distances of 8.14 Å. The determination of 24 couplings for a larger supercell as listed in Supplementary Table 3 fully confirms the three largest couplings, shown as empty symbols in Fig. 6b. It also shows that Y-kapellasite is a very two-dimensional material, and in the following, we neglect all interlayer couplings. Among the three second and six third nearest neighbor couplings of the distorted kagome lattice, the largest is J₇ with a value of 2% of J, which is rather small. This means that it is justified to focus the study of Y-kapellasite on the nearest-neighbor Hamiltonian. The couplings for the S_ND structure of Y-kapellasite are given in Supplementary Fig. 8 and Supplementary Table 2, respectively. There is one clear difference between S_XRD and S_ND: J_⬡ is 13% smaller than J for the S_XRD structure but 3% larger for the S_ND structure. We will discuss the implications for the Hamiltonian in the next section. We emphasize that, according to the ab initio calculations above, the value of the exchange coupling $J^{\prime}$ is considerably smaller than those of J_⬡ and J, which are of comparable size. In conclusion, by calculating and inspecting a large number of exchange couplings, we verified that the J, J_⬡, ${J}^{\prime}$ Hamiltonian is not a simplifying assumption but a defining feature of the material Y-kapellasite.

Fig. 6: Structure and exchange couplings of Y₃Cu₉(OH)₁₉Cl₈.

a Crystal structure of Y-kapellasite²⁶ (space group 148, $R\bar{3}$) with DFT relaxed hydrogen positions. b Exchange couplings of Y₃Cu₉(OH)₁₉Cl₈ as a function of interaction strength U, determined by energy mapping using GGA+U at J_H = 1 eV . Solid symbols: P1 cell, 7 couplings extracted. Empty symbols: $\sqrt{2}\times \sqrt{2}\times 1$ supercell, 24 couplings extracted. The inset shows the nearest-neighbor exchange paths of the perfect kagome lattice which differentiate into J_⬡, J, and $J^{\prime}$ in Y-kapellasite.

In what follows, we concentrate on the magnetic properties of Y-kapellasite as described by the structure S_XRD and the coupling constants J, J_⬡, $J^{\prime}$ of Supplementary Table 1 which place the material in the $\overrightarrow{Q}=(1/3,1/3)$ ordered regime of the classical phase diagram (see Fig. 2). This placement puts our calculations in agreement with the very recent experimental observation of (partial) magnetic order in Y-kapellasite³⁸.

Classical Monte Carlo simulations

We start our investigation of the Heisenberg Hamiltonian of Y-kapellasite (for the S_XRD structure) using the classical Monte-Carlo technique. Despite neglecting quantum fluctuations, this analysis allows us to study how thermal fluctuations impact the $\overrightarrow{Q}=(1/3,1/3)$ order (see Methods for technical details). In Fig. 5, we present results on the spin susceptibility in momentum space [see Eq. (7)] for different temperatures. At high temperatures [Fig. 5e], the response is almost homogeneous along the edges of the extended Brillouin zone. This response resembles the one of the standard undistorted nearest-neighbor kagome models, indicating that at these high temperatures details of the precise detuning between J, J_⬡, and $J^{\prime}$ do not yet affect the susceptibility. When T is lowered [going from panel e to panel a in Fig. 5], additional features become discernible such as three maxima around each corner of the extended Brillouin zone ($\sqrt{3}\times \sqrt{3}$ positions). Each such triad forms an equilateral triangle and with decreasing temperature the peaks become sharper. Simultaneously, the triangles show a slight rotation around their center points ($\sqrt{3}\times \sqrt{3}$ positions) until in the low-temperature limit the peaks reach the $\overrightarrow{Q}=(1/3,1/3)$ order positions [${K}_{2}^{\prime}$ points in Fig. 1b]. This shift of peaks roughly occurs along a line connecting the ${K}_{2}^{\prime}$ and Γ″ points. Please see Fig. 8 c for a trace of the peak positions at different temperatures. Note that as a result of the Mermin–Wagner theorem, real long-range magnetic order is possible only at strictly T = 0. However, the fact that the short-range correlations in the intermediate temperature regime manifest themselves in susceptibility peaks at incommensurate wave vectors away from $\overrightarrow{Q}=(1/3,1/3)$ indicates that thermal fluctuations act in a non-trivial and unexpected way. At T = 0, where the $\overrightarrow{Q}=(1/3,1/3)$ sets in, the maxima of the susceptibility are located at the ${K}_{2}^{\prime}$ points, with smaller peaks appearing at the $K_{1}^{\prime}$ and $K_{3}^{\prime}$ points. The relative heights of the latter peaks are ${\chi }_{K_{1}^{\prime} }/{\chi }_{{K}_{2}^{\prime}}\approx 40 \%$ and ${\chi }_{K_{3}^{\prime} }/{\chi }_{{K}_{2}^{\prime}}\approx 20 \%$, respectively.

Variational Monte Carlo results

We now analyze the ground-state properties of the spin model in the quantum regime with VMC calculations. As in the previous section, we focus on the set of exchange couplings obtained for the S_XRD structure of Y-kapellasite, which lies in the $\overrightarrow{Q}=(1/3,1/3)$ ordered region of the classical phase diagram (see Fig. 2). Our variational method is based on Gutzwiller-projected fermionic states (see Methods and Supplementary Note 5 for details). Optimizing the variational state, we obtain a finite value of the magnetic field variational parameter, which indicates the resilience of the classical $\overrightarrow{Q}=(1/3,1/3)$ magnetic order against quantum fluctuations. To corroborate this finding, we compute the spin susceptibility, Eq. (7), with 〈⋯〉 = 〈Ψ₀∣⋯∣Ψ₀〉 representing the expectation value over the optimal variational wave function. The results for a finite cluster of N = 972 sites (type-II, L = 6, cf. Supplementary Note 1) are shown in Fig. 7: the susceptibility is clearly dominated by sharp Bragg peaks at the ${K}_{2}^{\prime}$ points of the extended Brillouin zone, thus confirming the presence of $\overrightarrow{Q}=(1/3,1/3)$ magnetic order. We note that ${\chi }_{\vec{k}}$ is not significantly different from the classical result at zero temperature (i.e., no specific features are detected except for the Bragg peaks), despite the important contributions of the fermionic hoppings and the Jastrow factor to the variational energy. An almost identical susceptibility is obtained when considering the exchange couplings of the S_ND structure. Thus, according to our VMC results, the minimal Heisenberg model for Y-kapellasite has a $\overrightarrow{Q}=(1/3,1/3)$ magnetically ordered ground state.

Fig. 7: Variational Monte Carlo result.

Spin susceptibility in momentum space [Eq. (7)] computed by VMC. The results refer to the optimal variational state for the Hamiltonian of Eq. (1) with the exchange couplings of the S_XRD structure of Y-kapellasite (see Supplementary Table 1). The calculation is performed on a finite cluster of N = 972 sites (type-II, L = 6, cf. Supplementary Note 1).

Pseudofermion functional renormalization group calculations

Next, we employ the PFFRG approach^{40,41,42,43,44} to investigate ground state quantum effects in our distorted kagome Heisenberg model from a complementary methodological perspective. Within PFFRG, we compute the static spin susceptibility in momentum space ${\chi }_{\vec{k}}^{{{\Lambda }}}$ as a function of the RG parameter Λ (which acts as a low-energy frequency cutoff). We employ two variants of this technique, the one-loop and two-loop schemes, where the latter can be considered more accurate (but computationally more demanding) as it includes additional diagrammatic contributions to better account for the system’s fluctuations beyond mean-field (see Methods for details). Most importantly, an onset of magnetic order can be observed as an instability during the RG flow of the maximal $\vec{k}$-space component of ${\chi }_{\vec{k}}^{{{\Lambda }}}$. Such instability is indeed evident in the RG flows for both schemes (one-loop, two-loop) and for both structures of Y-kapellasite (see Fig. 8 a) confirming the findings from VMC. However, the fact that these instability features are quite weak and only detectable as small kinks rather than sharp peaks indicates the significance of quantum fluctuations, possibly associated with a small ordered moment. Note that the instability is observed at much smaller Λ in the two-loop scheme as compared to the one-loop approach, which is a known property resulting from the better fulfillment of the Mermin-Wagner theorem in the former method⁴¹. The momentum resolved susceptibility ${\chi }_{\vec{k}}^{{{{\Lambda }}}_{{{\scriptsize\mbox{crit}}}}}$ at the critical RG scale Λ_crit from two-loop PFFRG for the S_XRD structure is shown in Fig. 8b. The maxima are rather broad, again indicating strong effects of quantum fluctuations. Furthermore, the peaks do not exactly coincide with the $\overrightarrow{Q}=(1/3,1/3)$ positions, as in VMC results, but show a small shift along the ${K}_{2}^{\prime}-{{\Gamma }}^{\prime\prime}$-line, resembling our above findings from classical Monte Carlo. This indicates that quantum fluctuations may have similar effects as thermal fluctuations. We emphasize that this behavior is rather unusual, since, typically, quantum fluctuations tend to lock magnetic orders at commensurate positions. In Fig. 8c, we summarize the peak positions from one-loop and two-loop PFFRG as well as from classical Monte Carlo at intermediate temperatures. As can be seen, all results show a shift along similar momentum-space directions, however, the displacement away from the $\overrightarrow{Q}=(1/3,1/3)$ point becomes smaller as we advance the approach from one-loop to two-loop. It is hence conceivable that the shift would completely disappear upon further improving the method toward multi-loop schemes^45,46. We leave this as an open question for future investigations. We remark, however, that VMC and PFFRG both find magnetic long-range order for Y-kapellasite.

Fig. 8: Pseudofermion functional renormalization group results.

a One- and two-loop spin susceptibility flows from PFFRG for S_XRD and S_ND structures. The arrows mark a kink or peak during the flow indicating the onset of magnetic order. b Spin susceptibility in $\vec{k}$ space [Eq. (7)] from two-loop S = 1/2 PFFRG for S_XRD structure at the critical cutoff (marked by the red arrow in (a)). The extended Brillouin zone is indicated by the black hexagon. Maxima of ${\chi }_{\vec{k}}^{{{{\Lambda }}}_{{{\scriptsize\mbox{crit}}}}}$ appear at incommensurate positions. c Momentum-space position of the maximal susceptibility at the breakdown of the PFFRG flow together with the corresponding values obtained from classical Monte Carlo at finite temperature T. We show a section of the first Brillouin zone (gray lines) close to one edge of the extended Brillouin zone (black line), confer Fig. 1b. The PFFRG results for the S_XRD parameters are shown as triangles (the peak positions for the S_ND structure are similar). Red and green triangles represent one-loop and two-loop results, respectively. Dots represent the susceptibility maxima from classical Monte Carlo for the system with S_XRD couplings and the color represents temperature.

Linear spin-wave theory

We conclude the analysis of the magnetic properties of Y-kapellasite by showing in Fig. 3d the spin-wave spectrum and intensities for the Heisenberg Hamiltonian corresponding to the S_XRD structure. We observe that the spectrum is very similar to the simpler case of J = J_⬡ and $J^{\prime} =0$ (Fig. 3), which can then be regarded as a reliable minimal approximation for the full model. The intensity is largest at the ${K}_{2}^{\prime}$ point [see Fig. 1b] corresponding to the $\overrightarrow{Q}=(1/3,1/3)$ order, as also observed in the spin susceptibility results above.

Discussion

Summarizing, by a combination of DFT, effective spin models, classical (iterative minimization, classical Monte Carlo), and quantum approaches (VMC, PFFRG) we investigated the magnetic properties of a distorted kagome lattice as realized in the recently synthesized Y-kapellasite. We found an unexpectedly rich phase diagram already at the classical level which includes a collinear $\overrightarrow{Q}=0$ magnetic phase, two unusual non-collinear coplanar $\overrightarrow{Q}=(1/3,1/3)$ magnetic phases, and a classical spin liquid phase that resembles the jammed spin liquid phase found in the context of a bond-disordered kagome antiferromagnet. Our analysis of the spin model for Y-kapellasite places this system in the region of $\overrightarrow{Q}=(1/3,1/3)$ magnetic order with an excitation spectrum that lies halfway between that of an underlying triangular lattice of hexagons and a kagome lattice of trimers.

While it is not experimentally settled whether Y-kapellasite orders magnetically, our theoretical results provide strong evidence in favor of a magnetic $\overrightarrow{Q}=(1/3,1/3)$ ground state. The presence of an extended classical spin liquid phase in the vicinity of our DFT couplings sheds additional interesting light on this compound. Possibly, through external perturbations such as pressure or strain one might be able to shift the couplings towards the classical spin liquid phase, which, due to the large extent of this regime, may not require any fine-tuning. This opens the question about the fate of the classical spin liquid upon adding quantum fluctuations, which we did not tackle in this work. Given the complexity of this phase already on the classical level, one may expect even richer phenomena in the quantum case, including a quantum spin liquid. The numerical investigation of this regime in the quantum limit will certainly be a challenging future task but also gives hope for rewarding insights. In total, this work demonstrates that a relatively simple but realistic distortion of the kagome lattice gives rise to a multitude of interesting and unexpected magnetic phenomena whose full investigation goes far beyond the scope of the present work. In the future, our investigation may inspire and guide both a deeper experimentally motivated investigation of Y-kapellasite, as well as a closer numerical analysis of the underlying spin model.

Methods

DFT-based energy mapping

We calculate the electronic structure of Y-kapellasite with DFT using the full potential local orbital basis set⁴⁷ and the generalized gradient approximation (GGA) to the exchange-correlation functional⁴⁸. We apply the GGA+U approximation⁴⁹ to correct for strong electronic correlations of the Cu 3d electrons. We set the Hund’s rule coupling to a typical value value^50,51J_H = 1 eV for Cu²⁺ and vary only the onsite interaction U. Even though the S_ND structure³⁷ nominally has 8/9 filling, there is no evidence that Y-kapellasite is charge doped, and therefore we treat the O1 position as occupied with a hydroxy group (or a chloride ion which leads to the same results). In this position, the S_XRD structure²⁶ has an orientationally disordered OH⁻ ion between two Y³⁺ ions, and therefore a 1/6 occupation of the six symmetry equivalent H positions is consistent with the $R\bar{3}$ space group. We model the orientationally disordered OH⁻ ion using the virtual crystal approximation⁵², setting the nuclear charge of H in this position to 1/6. The hydrogen positions H2 to H4 are relaxed within GGA in both structures. We shift the partially occupied H1 hydrogen position to the equilibrium O–H distance. The resulting structure is shown in Fig. 6a.

We use total energy mapping^53,54 to determine the Heisenberg Hamiltonian parameters of Y-kapellasite. For that, we calculate with DFT(GGA+U) the total energy for 24 out of the 47 unique spin configurations which are possible with the 9 inequivalent Cu²⁺ ions in the P 1 unit cell of Y₃Cu₉(OH)₁₉Cl₈. Considering that third-neighbor couplings are important for some kapellasite type compounds⁵⁰, we also perform calculations for a $\sqrt{2}\times \sqrt{2}\times 1$ supercell with 18 independent Cu sites; we calculate 44 out of nearly 30,000 spin configurations with distinct energies.

Iterative minimization

To numerically determine the classical ground state of the spin model of Eq. (1), we employ the iterative minimization method^55,56. We initialize our system in a random configuration and we iteratively perform local moves to update the spins. In each move, we pick up a random spin, ${\overrightarrow{S}}_{i}$, and we orient it antiparallel to the local field created by the neighboring spins, i.e.,

$${\overrightarrow{S}}_{i}\mapsto -\frac{{\overrightarrow{B}}_{i}}{| {\overrightarrow{B}}_{i}| },\quad \,{{\mbox{with}}}\,{\overrightarrow{B}}_{i}=\mathop{\sum}\limits_{j}{J}_{ij}{\overrightarrow{S}}_{j}.$$

(8)

The procedure is repeated for a sufficient number of steps until the energy converges. In order to reduce the risk of ending up in local minima, we repeat the calculations several times starting from different spin configurations and we keep the solution with the best energy. The calculations are performed on the small finite-size clusters shown in Supplementary Fig. 1 with N = 27 and N = 36 sites, and periodic boundary conditions. It is important to emphasize that finding a classical ground state with ${\overrightarrow{L}}_{\bigtriangleup }=0$ on one of these small clusters (with periodic boundary conditions) implies that one can immediately define a ${\overrightarrow{L}}_{\bigtriangleup }=0$ ground state for any larger cluster of the same type (see Supplementary Note 1).

Classical Monte Carlo simulations

We perform a Monte Carlo analysis using the Metropolis algorithm with the over-relaxation protocol for better thermal convergence^57,58,59,60. For the investigations of the S_XRD structure, the system that we simulate is a cluster of N = 7803 spins with periodic boundary conditions (type-II cluster with L = 17 (Supplementary Fig. 1). It is seeded at ${T}_{0}=2\ {J}_{\max }$ with random spins and cooled down via T = T₀e^−γn where n = 0, 1, 2, … is the number of Metropolis steps. During each step, every spin is updated once on average by a new random spin. The update takes place either with certainty if the acquired energy ΔE ≤ 0 or with a probability p = e^−ΔE/T. For the different coupling regimes, we set γ = 0.001 and cooled 100 (10) random systems down to $T=0.1\ {J}_{\max }$ ($T=0.001\ {J}_{\max }$), where ${J}_{\max }=\max ({J}_{\hexagon},J,J^{\prime} )$. For the investigation of the classical phase diagram, we restrict ourselves to type-II clusters with L = 5 (N = 675 spins), for a numerical speedup.

Linear spin-wave theory

We performed linear spin-wave calculations with the SPINW package⁶¹, computing the classical ground state by energy minimization.

Variational Monte Carlo

We employ Gutzwiller-projected fermionic states as variational ansätze⁶². This class of wave functions has been shown to provide an accurate description of the ground state of several spin models⁶³, including state-of-the-art results for kagome lattice antiferromagnets^{10,29,64,65,66}. The optimal variational ansatz for the Hamiltonian of Y-kapellasite is a Gutzwiller-projected Jastrow–Slater wave function possessing $\overrightarrow{Q}=(1/3,1/3)$ magnetic order (see Supplementary Note 5 for the definition). For the optimization of the variational parameters, we use the stochastic reconfiguration method⁶⁷.

Pseudofermion functional renormalization group

The PFFRG method is based on the one-loop plus Katanin truncation PFFRG scheme first introduced in ref. ⁴⁰ and extended to the two-loop plus Katanin variant in ref. ⁴¹. It utilizes the Abrikosov pseudofermion representation of spin operators. This spin representation enlarges the Hilbert space by adding two unphysical S = 0 states per site (unoccupied, doubly occupied) which, however, leave the ground state properties largely unaffected⁴². Within PFFRG the bare propagator of the fermions is regularized by a sharp cutoff function

$${G}_{0}(\omega )=\frac{1}{i\omega }\quad \longrightarrow \quad {G}_{0}^{{{\Lambda }}}(\omega )=\frac{\theta (| \omega | -{{\Lambda }})}{i\omega }\ .$$

(9)

Here, ω is a continuous Matsubara frequency at T = 0 and the cutoff Λ prohibits fermionic propagation if ∣ω∣ ≤ Λ. This insertion causes a cutoff dependence of the generating functional for the fermionic one-particle-irreducible vertex functions. Flow equations that describe the Λ derivatives of all n-particle vertex functions can be derived. These equations couple the n-particle vertex to all m-particle vertices with m ≤ n + 1 leading to an infinite hierarchy of equations. In principle, physical results in the cutoff-free limit Λ → 0 can be obtained by solving the integrodifferential flow equations starting from the limit Λ → ∞ where the initial conditions are set by the bare interactions from our spin Hamiltonian. For numerical solvability, this hierarchy of equations needs to be truncated. In the one-loop scheme, the truncation occurs on the level of the three-particle vertex which is replaced by contributions from the Katanin scheme⁶⁸, particularly, the single-scale propagator is upgraded to ${S}^{{{\Lambda }}}(\omega )=-\frac{d}{d{{\Lambda }}}{G}^{{{\Lambda }}}(\omega )$ where the full Green’s function is ${G}^{{{\Lambda }}}(\omega )={\left[{\left({G}_{0}^{{{\Lambda }}}(\omega )\right)}^{-1}-{{{\Sigma }}}^{{{\Lambda }}}(\omega )\right]}^{-1}$. The one-loop flow equations for the self-energy Σ^Λ and the two-particle vertex Γ^Λ are depicted diagrammatically in Supplementary Fig. 7. In the two-loop approach, further contributions of the three-particle vertex are included, which have the form of nested one-loop diagrams⁴¹. We solve the flow equations numerically with an Euler scheme in real space taking into account finite spin correlations on hexagonal clusters with an edge length of N ≥ 7 nearest-neighbor distances around reference sites from each sublattice. The Matsubara frequencies are discretized using a linear plus logarithmic mesh with M_ω ≥ 60 points. We carefully analyzed that the qualitative PFFRG results are converged with respect to the number of frequency points and the finite correlation length. From the resulting two-particle vertex, we are able to compute the Λ dependent static spin susceptibility ${\chi }_{\vec{k}}^{{{\Lambda }}}$ in momentum space. For more details, we refer the reader to refs. ^40,41,42.