Spiral-spin-liquid behaviors and persistent reciprocal kagom\'{e} structure in frustrated van der Waals magnets and beyond

We study classical $J_1$-$J_2$ models with distinct spin degrees of freedom on a honeycomb lattice. For the XY and Heisenberg spins, the system develops a spiral spin liquid (SSL) that is a thermal cooperative paramagnetic regime with spins fluctuating around the spiral contours in the momentum space, and at low temperatures supports a vector spin-chirality order despite the absence of long-range magnetic order. In a strong contrast, for the Ising moments, the low-temperature spin correlation forms a reciprocal"kagom\'e"structure in the momentum space that resembles the SSL behaviors and persists for a range of exchange couplings. The unexpected emergence and persistence of the reciprocal"kagom\'e"are attributed to the stiffness of the Ising moments and the frustration. At higher temperatures when the thermal fluctuations is strong and the spin correlation is not fully melted, the reciprocal structures evolve from the"kagom\'e"towards the ones demanded by the soft spin limit. This contrasts strongly with the behaviors of the spiral contours in the SSL regime for the continuous spins. We suggest various experimentally relevant systems including van der Waals magnets such as the transition metal phosphorus trichalcogenides TMPX$_3$, Cr$_2$Ge$_2$Te$_6$, the rare-earth chalcohalides (like HoOF, ErOF and DyOF) and other isostructural systems to realize the SSL-like behaviors and/or the reciprocal kagom\'{e} structure.


I. INTRODUCTION
The recently developed van der Waals (vdW) materials provide an excellent platform for the understanding of the two-dimensional physics and the potential application of various devices [1][2][3][4][5]. The vdW materials are three dimensional but due to the weak van-der-Waals force between the adjacent layers, a monolayer of vdW materials can be obtained through various exfoliation methods [5]. The transition-metal phosphorus trichalcogenides, TMPX 3 , are a class of vdW materials where the transition metals (TM) are combined with phosphorus (P), and chalcogenides (X = S, Se, Te). In such materials, TMs constitute a honeycomb lattice with intrinsic magnetism, and most members of the family exhibit an antiferromagnetic exchange [6]. Apart from their excellent structure, the magnetic vdW materials own their advantage in terms of their variety and the controllability. The spin Hamiltonian differs from material to material. The spin degrees of freedom in these materials could be of the Ising, XY or Heisenberg types [1,7,8]. For some materials, e.g. Cr 2 Ge 2 Te 6 , the spin type can even be tuned by applying the hydrostatic pressure [9][10][11]. Moreover, it is convenient to vary the anisotropic interactions by the external perturbations, such as the gating and strain, or the proximity effects [1,11]. As far as we are aware, in application, most of the current efforts are devoted to the realization of the long-range magnetic orders on the vdW materials as the magnetic orders are partially forbidden by the Mermin-Wagner theorem, and this may be used in * gangchen@hku.hk designing magnetic devices. Other efforts have been devoted to exploring the interesting magnetic excitations, such as the topological magnons [12,13] for the honeycomb lattice antiferromagnet CrI 3 , with respect to the magnetically ordered ground states [14,15].
Here, we deviate from the practical purpose of device designing with vdW magnets, and instead address the possibility of interesting fundamental physics that could potentially occur in these new materials. The direction that we are toward here is magnetic frustration. Frustrated magnetism has attracted tremendous interest for decades because of many unconventional and exotic properties [16,17] and the potential application to quantum computing and quantum information. Under thermal or quantum fluctuations or both, exotic states, i.e. classical or quantum spin liquids [16], spin ice [18][19][20], Kitaev spin liquid [21], Berezinskii-Kosterlitz-Thouless (BKT) phase [22] and many others could emerge. Frustration can come from the lattices themselves. The well-known lattices of geometric frustration are the triangular lattice, the kagomé lattice, and the pyrochlore lattice. Generally speaking, frustrations come from the competing interactions, e.g. the competition between the inter-sublattice and the intra-sublattice interactions. The antiferromagnetic J 1 -J 2 spin model on a square lattice is a simple example of competing interactions [23]. Another wellknown example is the same model but on a diamond lattice [24]. In this model, there exists a spiral spin liquid (SSL) regime within some parameter region and this physics is already detected in experiments in the diamond lattice antiferromagnet MnSc 2 S 4 [25]. The SSL is a special family of classical spin liquids [18,20,24,[26][27][28][29][30][31][32][33][34], and it exhibits partial degeneracies where, in the thermodynamic limit, the spin structure factor displays the spiral surfaces or spiral contours in the reciprocal space [24,[35][36][37][38][39][40][41][42][43][44]. To our knowledge, there are relatively limited numbers of works about the SSL physics, and most previous studies have focused on the Heisenberg spins. The twodimensional (2D) vdW magnets provide an opportunity to explore the physics of the SSLs and render new magnetic degrees of freedom beyond the Heisenberg spins for the study of the SSLs and other frustrated spin physics.
Apart from these transition metal-based-vdW magnets, there has been intensive interest in the rare-earth magnets. Recently, the rare-earth honeycomb lattice magnets have been proposed as candidate for Kitaev materials [45][46][47], and a series of vdW rare-earth chalcohalides with an equivalent honeycomb geometry have been synthesized [48]. Anisotropic interactions are quite common for the rare-earth magnets [47], and one such anisotropic limit is the Ising model when the spin moment is Ising-like. The Heisenberg model is applicable to the Gd-based magnet with S = 7/2 and was also argued to be relevant for some Yb-based magnets (where the moment is effective spin-1/2) [49]. Since many 2D vdW magnets have a 2D honeycomb structure, we consider a J 1 -J 2 spin model on a honeycomb lattice with the Ising, XY or Heisenberg spins. For the XY and Heisenberg spins, Refs. [24,36,39,42] have shown that there exist SSLs on bipartite lattices such as the honeycomb structure. The magnetic-order transition, which breaks the U(1) or SO(3) symmetry, will not happen in the 2D systems according to the Mermin-Wagner theorem. Nevertheless, we find that the Z 2 symmetry, i.e. the chiral symmetry, is spontaneously broken at low temperatures when J 2 is larger than a critical value ∼ 0.7J 1 .
Although the continuous spin seems necessary for the construction of the spin spirals, a pattern similar to the spiral contour for the SSL, a reciprocal "kagomé" structure, will emerge for the low-temperature spin correlation in the momentum space for the Ising spin. This behavior resembles the SSL for the continuous spin. However differing from the varying spiral contour structures with varying parameters, our calculation in Sec. III shows that this reciprocal "kagomé" structure persists at low temperatures for a range of J 2 's. This remarkable result, as we further explain in Sec. III, is a unique but rather natural property of the Ising spin moment and arises from the stiffness of the Ising moment and the local constraint due to the frustration upon varying of the parameters. It is also found interesting to explore the thermal evolution of the reciprocal structure. Unlike the robustness against the variation of J 2 , the reciprocal structure is found to deviate from the reciprocal "kagomé" as the temperature is above a crossover temperature in Sec. IV. The reciprocal structure gradually evolves toward the contours that are demanded by the exchange interaction with the given exchange couplings in the soft spin limit. This thermal behavior at high temperatures is further understood by the soft spins due to the thermal fluctuations.
The rest of the paper is organized as follows. In Sec. II, we explain the model and the basic properties of the SSLs (a) for the continuous spins. Some of the physics is explained from the local-constraint point of view. While the case of continuous spin has been studied in previous work, the perspective of the local energetic constraint provides some insights into the emergent properties. Moreover, we show the appearance of the finite-temperature spin chirality order that relates to the local electric polarization via the inverse Dzyaloshinskii-Moriya mechanism. In Sec. III, we turn our attention to the Ising spin where the SSL-like behaviors surprisingly occur and persist for a range of parameters at low temperatures. Apart from the finite-temperature phase transition due to the discrete nature of the local moments, the spin correlation supports the reciprocal "kagomé" structure in the momentum space at the low temperatures. This is explained from the local constraint point of view. In Sec. IV, we further explore the temperature evolution of the spin correlation for the Ising spin and explain the peculiar crossover behaviors in the reciprocal structures. A comparison with the self-consistent Gaussian approximation is also given. Finally in Sec. V, we summarize our re- sults and discuss the experimental relevance. We particularly emphasize the Ising spin connection with the rareearth chalcohalides (such as HoOF, ErOF and DyOF) and make predictions based on our theoretical results.

II. MODEL WITH CONTINUOUS SPINS
We consider the J 1 -J 2 spin model on the honeycomb lattice with periodic boundary conditions ( Fig. 1), where ij ([ij]) refers to the nearest (next-nearest) neighbors and S i is the classical Ising, XY or Heisenberg spin on site i. The exchange coupling J 2 is antiferromagnetic, and J 1 could be ferromagnetic or antiferromagnetic because a simple transformation on one sublattice switches the sign of J 1 . Here we set J 1 to be antiferromagnetic. With only J 1 , the ground state is a simple Néel state because of the bipartite nature of the honeycomb lattice. With only J 2 exchange, the Hamiltonian decouples to two independently antiferromagnetic spin models on two triangular sublattices. For the Ising spins, the ground state is massively degenerate, and it is 120 • state for the XY and Heisenberg spins. When J 1 and J 2 are both non-zero, a large frustration is introduced when J 2 is larger than a critical value J 2c and interesting properties could appear due to the strong frustration. To obtain comprehensive and accurate behaviors of this model, we mainly employ the classical Monte Carlo (MC) method to investigate both the zero-and finite-temperature properties for three different types of the spin moments.

A. Analytical results for continuous spins
In this section, we focus on the XY and Heisenberg spins for which the model has a global U(1) and SO(3) symmetry, respectively. Thus, there is no finitetemperature magnetic ordering transition according to the Mermin-Wagner theorem. Although part of the results in this section was previously known in Ref. [36], which studied the quantum Heisenberg spins, we include this analysis for completeness and for later comparison with the peculiar Ising case. In the Luttinger-Tisza method, the local constraint |S i | 2 = S 2 for each spin is softened to a weak global constraint i |S i | 2 = N S 2 , where N is the number of lattice sites. In this weak constraint, we minimize the energy and check whether the strong constraints are satisfied afterwards. If these strong constraints are satisfied as well, the ground state obtained from the Luttinger-Tisza method is just the ground state of the initial model. In practice, we define the Fourier transformation of S i in the sublattice µ = A, B as S µ i = (N/2) −1/2 k e ik·ri S µ (k). The Hamiltonian in Eq. (1) can be rewritten as with φ(k) = (S A (k), S B (k)) T . Here J (k) is the Fourier transform of the adjacency matrix on the honeycomb lattice and is given as where the vectors b µ refer to the nearest neighbor vectors of the honeycomb lattice. The minimal eigenvalue of matrix J (k) is given as with and the energy of the ground state is the minimum value of − (k)/2. The range of Λ(k) is from 0 to 3, and γ c = (J 2 /J 1 ) c = 1/6 is a critical point. When γ ≡ J 2 /J 1 < γ c , the point of Λ(q) = 3 will minimize − (q), i.e., q = 0, and this corresponds to the Néel order.
For larger values of γ, − (q) takes the minimum when Λ(q) = (2γ) −1 , and these momentum vectors constitute closed contours in the reciprocal space, dubbed the spiral contours, as shown in Fig. 2. Clearly, the system has a massively degenerate ground-state manifold, indicating a strong frustration. With increasing γ, the spiral contour expands around the Γ point and touches the M point when γ = 1/2. For a larger γ > 1/2, a single spiral contour splits into several contours around the K points in the first Brillouin zone and they gradually shrink to K points. At the limit γ → ∞ or J 1 = 0, the spiral contours disappear leaving only a single spiral state, the 120 • state. The model reduces to an antiferromagnetic model on the triangular lattice. An intuitive method of obtaining the information about the ground state is a geometric one. This method is often adopted in frustrated systems to find a local constraint to minimize the classical (and sometimes even quantum) ground-state energy. In the classical spin ice on a pyrochlore lattice, for example, the ice rule is a local constraint that should be satisfied in the ground-state manifold [26]. The key ingredient of this trick is to split the Hamiltonian into equivalent cluster units, and then one minimizes the energy of each cluster unit to obtain the lowest energy and acquire the local constraints for each cluster unit. We rewrite Eq. (1) as where S = S 1 + S 2 + S 3 is the sum of three corner spins on a unit / and S 0 is the central spin on the unit. In the following, we will use to represent both and . Figure 1(a) shows the honeycomb lattice and a unit. In a honeycomb lattice with N sites, there are N 's. From Eq. (6) one sees that γ c = 1/6 is a critical value and is the same as the Luttinger-Tisza result because the largest length of S is three times the length of S 0 . If γ < γ c , the minimal energy is reached when S = −3S 0 , i.e. S 1 , S 2 , and S 3 are all antiparallel to S 0 . This state is simply the Néel state. In the region of γ > γ c , the condition, for each unit can always be satisfied for the XY and Heisenberg spins to minimize the energy. In Fig. 1(b), we depict the evolution of the spins on with increasing γ. We find that there exists a degeneracy of each unit under the local constraint S 0 + 2γS = 0 when γ ≥ γ c . This leads to the massive degeneracy of the ground state. In the limit γ = ∞ or J 1 = 0, the constraint on each unit reduces to the weak constraint of S = 0 for each and the massive degeneracy is lifted to the discrete Z 6 degeneracy, i.e., the ground-state order is the 120 • state. In addition to the continuous rotational symmetry breaking, this state also breaks the Z 2 symmetry or the chiral symmetry that describes the spin rotation pattern on a triangle clockwise or counter-clockwise. This symmetry is discrete and can be spontaneously broken at finite temperatures, which means that there can be a finitetemperature chiral transition under this limit [50,51]. Moreover, we expect that as long as γ is large enough, this discrete symmetry can still be spontaneously broken at a finite temperature and the chiral order will occur.
The other important quantity is the spin structure factor and it can be detected experimentally [24]. The lowtemperature spin structure factor provides an important characterization of the physical properties related to the classical ground state degenerate manifold. In this paper, we define the spin structure factor on the A or B sublattice as with µ = A or B. In zero temperature, only with those q's minimized Eq. (4), S µµ k are nonzero. The corresponding spin structure factors can be expressed as In Fig. 2, we plot the spin structure factors for different γ's where the degenerate momentum vectors form the spiral contours. When 1/6 < γ < 1/2, the spiral contour is a single closed loop around the Γ point in the first Brillouin zone, and its size becomes larger with increasing γ. This contour touches the first Brillouin zone boundary at the M point when γ reaches 1/2. It splits to several contours around the K points when γ > 1/2. In the limit of γ → ∞ or J 1 = 0, these spiral contours shrink to the K points, which indicates the rise of the 120 • state.

B. Numerical simulation and finite temperature thermodynamics
Here, we use a Metropolis algorithm combined with the over-relaxation method [52,53] and the parallel tempering [54] to simulate the proposed model at finite temperatures. As we are all aware, frustrated systems generally have a large energy barrier between numerous local minimal energy states with only local updates such as the Metropolis updates. This will lead to spin configurations that fluctuate near the minima for a long time and make the simulation unreliable. To overcome this obstacle, a simple and efficient approach is to use the parallel tempering scheme.
In the parallel tempering scheme, multiple replicas of the same system, randomly initialized, are simulated at different temperatures. Then the exchange of replicas between the nearest temperatures occurs with a certain probability, and the exchange-acceptance ratio is calculated according to the detailed balance condition. Replicas can swap around the whole temperature region and this significantly suppresses the configuration freezing for replicas at low temperatures. This is because updates are efficient for high-temperature replicas and these spin configurations can gradually convey to the replicas of the low temperatures. The key to this approach is to select appropriate temperatures for every replica to ensure that the exchange-acceptance probabilities between the nearest replicas are not too small. In this paper, we carefully select the temperatures for each replicas to make the exchange probabilities not smaller than 0.32. For this goal, we first follow the feedback-optimized plan in Ref. [54] to obtain the tentative temperatures and the corresponding energies for a small system size. According to these temperatures and energies, temperatures of replicas for other system sizes can be calculated with a given exchanged-acceptance ratio. The reason that this strategy works is that the energy density is almost independent of the system sizes and the energy-temperature curve is continuous in this model. Moreover, we employ the over-relaxation method, which could improve the performance of the simulations for the continuous spins. We carry out 128 independent simulations, and each one contains multiple replicas of the same system but in different temperatures. In every independent simulation, a whole MC step consists of a single Metropolisupdate sweep and subsequently a single over-relaxation sweep. Besides, a parallel-tempering update will occur after every 50 MC steps. After thermalizing systems to equilibration, 4 × 10 4 samples are produced in each sim-ulation, and in total 5.12 × 10 6 samples are used for the statistical analysis.
To determine the finite-temperature phase diagram, we measure the specific heat C v and the chiral order parameter m c . The peaks of the specific heats C v can be used to determine the crossovers or the phase transition. The chiral order parameter m c on a sublattice lattice is defined as where t is taken over all up-triangular units in the sublattice A or B and ij∈t sums over three bonds of each up triangle t in a clockwise order. In Figs. 3(a-c), we plot the curves of the specific heat C v and the chiral order parameter m c at γ = 0.08, 0.40, 0.80 with the linear system sizes L = 12, 24, 36 for the XY spin. When γ = 0.08, the system undergoes a crossover from the quasi-Néel state of the low temperatures to the paramagnetic state of the high temperatures, and the specific heat only has a non-divergent round peak at the crossover due to the Mermin-Wagner theorem. After entering the SSL regime for γ = 0.40, there exists a crossover accompanied by a nondivergent round peak without chiral-symmetry breaking. For a larger nextnearest exchange interaction, γ = 0.80, the frustration is not large enough to prevent the appearance of the spontaneous breaking of the chiral symmetry. This leads to a sharp peak in the specific heat but the magnetic order is still forbidden by the thermal fluctuation. At the same time, the chiral order parameter m c has a rapid growth near the transition indicating the breaking of the chiral symmetry. Using the highest-temperature peaks of the Due to the complex behaviors of the specific heat near γ = 0.40, the highest-temperature peak is taken as the phase boundary or crossover. Comparing the phase boundaries or transitions of different system sizes, we find that L = 36 is large enough to determine the phase boundaries. The vertical line at γc = 1/6 separates the unfrustrated and frustrated regions. Near γ = 0.70, the thermal crossover switches to the phase transition for γ > 0.70 where the chiral symmetry is spontaneously broken. This approximate γ region is marked with a light orange band. specific heat C v and the points of the rapid increase in the chiral order parameter m c , the phase diagram of the XY spin can be determined, as shown in Fig. 4(a). Near γ ≈ 0.7, the crossover changes to a phase transition. In the Heisenberg case, the same procedure is applied as shown in Fig. 3(d-f), and Fig. 4(b) is the phase diagram. Due to the system size and the numerical method, the potential BKT transition is not discernible here. For both the XY and Heisenberg spins, the behaviors of specific heats at γ = 0.40 are strongly affected by the system sizes. We think this is due to the stronger frustration. Near γ = 0.25, the energy difference between the Néel order and the spiral spin states with the wavevectors from the spiral contour is small and at finite temperatures, there exists strong competition between these spin configurations. This leads to rather complex behaviors and makes the simulations more difficult.
One outcome of the vector spin chirality order at low temperatures is the local electric polarization. This is obtained from the inverse Dzyaloshinskii-Moriya mech-anism that gives the local electric polarization [55], P ij ∝ e ij × (S i × S j ), where e ij is the unit vector that connects site i to site j. A finite vector spin chirality order implies a distribution of the local electric polarization P ij that could lead to a modification in the electric response and the structural distortion.

C. Spin structure factors
Here we further determine the evolution of the spin structure factors at different γ's for the XY and Heisenberg spins. The spin structure factors can be detected by the neutron scattering to reveal the magnetic structures. We adopted three different γ's and three different temperatures for both spins. The results are presented in Fig. 5.
The temperatures in Fig. 5 are marked in Fig. 4 with the purple points except for T 1 = 3J 1 . We choose T 1 > T 2 > T 3 in the plots. In the following, we analyze the case of the XY spin and it is similar for the Heisenberg spin. In Fig. 5(a 1 -c 1 ) of the XY spin, we set γ = 0.08, and the peaks of the spin structure factors are always at the Γ point as expected for the proximate Néel state at zero temperature. At the high temperature T 1 , the system is in the disordered state with round peaks; the system is in the quasi-Néel state at the low temperature T 3 with a sharpened peak. In Fig. 5(d 1 -f 1 ), we plot the spin structure factors for γ = 0.40. At T 1 , there are broad peaks located at the K points which means the spin structure factors are dominated by the J 2 exchange. At temperature T 2 , the peaks of the spin structure factors form a visible spiral contour and there are massively degenerate states. Further decreasing the temperature to T 3 , the spiral contour becomes more sharp accompanied by some peaks that will be discussed next. For γ = 0.8, the spin structure factors of the high temperatures in Fig. 5(g 1 ) cannot be distinguished from those of Fig. 5(d 1 ). Nevertheless, there are several spiral contours around the K points for the low temperatures T 2 and T 3 , as shown in Fig. 5(e 1 ,f 1 ). At the temperature T 3 , the discrete spiral contours are more clear. For the Heisenberg spin, the numerical results of the spin structure factors are shown in Fig. 5(a 2 -i 2 ).
In Fig. 5(e 1 ,f 1 ,h 1 ), it seems that there are magnetic order peaks. According to the Mermin-Wagner theorem, however, this model should not have any magnetic order at finite temperatures. The thermal order-by-disorder mechanism to maximize the entropy work similarly as the three dimensions [24,42]. Near these wavevectors where the peaks are located, the entropies are larger than those of remaining parts on the spiral contours, but thery are not strong enough to induce any long-range magnetic order. With further decreasing temperature, this entropic effect becomes weaker and ultimately, at zero temperature, the relative difference disappears as shown in Fig. 2(b,c) with an equal strength. To give a consistency check of the absence of the true long-range order, we simulate the system at temperatures lower than the crossover temperature for γ = 0.5 where the peaks of the spin structure factor are located at the M point. The M point magnetic order on the triangular lattice was previously known to be a stripe-type magnetic structure where the spins are ferromagnetically aligned along one Bravais lattice vector direction and are antiferromagnetically aligned along the other Bravais lattice vector direction [56]. The results for the spin correlation and the spin structure factors are summarized in Fig. 6.
The spin correlation C(r) = S(r) · S(0) in Fig. 6 is measured along one Bravais lattice vector direction and the distance r is chosen to be an even lattice spacing. C(r) has no sign oscillation and can reflect whether there is magnetic order. In the XY case ( Fig. 6(a)), the spin correlation decays exponentially with increasing distance. For the Heisenberg spin, C(r) also decays as the distance r increases, as can be seen from Fig. 6(b). What is different from the XY spin is that the spin correlation C(r) is more like a power-law function. This discrepancy makes the peaks of the spin structure factor in the Heisenberg spin more pronounced than those in the XY spin as shown in Fig. 6. Towards large distances r, the deviation from the fitted line is a consequence of the pe- riodic boundary conditions and, with larger system sizes, more points will fall on the fitted line. We do not give much attention to the spin structure factors in the high-temperature paramagnetic state. In such a regime, the key magnetic property is submerged by the strong thermal fluctuations. At the infinite temperature limit, this is totally true but at not high temperature, thermal fluctuation may not erase all information. Comparing Fig. 5(a 1 ) with Fig. 5(d 1 )(g 1 ), they are obviously different, though they are all in the paramagnetic phase with the same temperature T 1 = 3J 1 . This difference originates from different values of γ. In the above context, we know that at zero temperature γ < γ c and γ > γ c drive the system into the Néel state and the spiral spin liquid, respectively. For γ = 0.08, the spin configurations of or near the Néel order have lower energies than other configurations. At high temperatures, the peaks of spin structure factors will broaden around the Γ point. For other γ values, the peaks of the spin structure factors for the ground state will broaden at high temperatures as well. Qualitatively, Fig. 5(d 1 ) and Fig. 5(g 1 ) are the results of spiral contours in Fig. 5(f 1 ) and Fig. 5(i 1 ) after thermal broadening, respectively. Although most information about the low-temperature property is already drowned in the thermal fluctuation, the difference between Fig. 5(a 1 ) and Fig. 5(d 1 )(g 1 ) still can tell us the spin structure factors rough shape at low temperatures. That is to say, if a similar structure, such as that in Fig. 5(a 1 ), is measured in experiment, it is enough to exclude the probability of the SSLs.

III. PERSISTENT RECIPROCAL "KAGOMÉ" STRUCTURE FOR ISING SPIN
For the Ising spin, one cannot construct the (noncollinear) spin spirals with the Ising moments. As we demonstrate below, however, there exists an emergent "kagomé" structure of the spin correlation in the reciprocal space due to the Ising nature of the local moment and the frustration. This reciprocal "kagomé" structure is reminiscent of the spiral contour at γ = 0.5 for the continuous spins. In addition, we show that this structure is persistent for the Ising spins as long as γ > 1/4. We attribute this phenomenon to the "stiffness" of the Ising spin.
To begin with, we first rely on the geometric method to obtain the lowest energy and the constraint of groundstate spin configurations for the Ising spin. Due to the discrete nature of the Ising spin, the upper limit of γ for the Néel order is no longer 1/6, which is the value for the continuous spins. The spin configuration of the Néel order is the same as the continuous spin such that S 1 , S 2 , and S 3 are all antiparallel to S 0 with only the global Z 2 symmetry (see the cluster unit in Fig. 1). If flipping one of S 1 , S 2 , or S 3 reduces the total energy compared with the Néel state, this indicates that, the Néel order is no longer the ground state, and the critical value of the Néel order is found to be γ c = 1/4. When γ > γ c , the ground state is massively degenerate as shown in Fig. 1, and the local constraint is two up spins and two down spins on each unit . Remarkably, we find that this constraint stays the same and valid for any γ > γ c , and this means that the spin structure factors of the ground state for different γ's share similar structures. This is fundamentally different from the case for the continuous spin where the spiral contour varies with the parameter γ. In the following, we perform the numerical simulation and provide the theoretical understanding. It is found that, the "stiffness" of the Ising spin pins the spin structure factors at low temperatures to a reciprocal "kagomé" structure in the momentum space. Moreover, the Z 2 symmetry will not lead to the chiral symmetry and the chiral phase transition does not exist for any γ.
Employing the same Monte Carlo simulations but without invoking the over-relaxation method, which is invalid for discrete spin, we show the specific heat C v at γ = 0.10, 0.40, 0.80 in Fig. 7. The finite-temperature phase diagram can be obtained as shown in Fig. 8 where the "boundary" is indicated by the peaks of the specific heat. When γ < γ c , the phase boundary is the Néelparamagnetic phase transition where the specific heat is divergent accompanied by spontaneous Z 2 symmetry breaking. In spite of there being no theorem to restricting the occurrence of the long-range orders, there is only crossover at finite temperatures for γ > γ c due to strong frustration. The sharp peaks in the specific heat occur at lower temperature for larger system sizes, and only round maxima are expected to occur at high temperatures for the thermodynamic limit with L → ∞. The phase diagram in Fig. 8 shows that the crossover temperatures increase with γ as at this time the more important interaction comes from the exchange J 2 . The round maximum shown with the dashed curves of Fig. 8 is a crossover separating the less correlated higher-temperature regime and more correlated low-temperature regime where the reciprocal kagomé regime becomes more and more visible. In highly frustrated magnets, a round maximum at the finite temperature of the specific heat is quite common and often observed. It is simply an indication of the large entropy loss at the crossover temperature point.
Previous works on the same Ising model also explored the thermodynamic properties. In Ref. 57, the authors gave a very similar phase diagram of the model via unsupervised machine learning with the system size fixed to 900 sites based on training data obtained from Monte Carlo simulations with the Metropolis update and thermal annealing. In addition, the authors of Refs. 58 and 59 did interesting work on the same model with some focus on the parameter region with γ < γ c where the ground state is a simple Néel state. Through finite-size scaling analysis, they found that the phase transition belongs to the 2D Ising universality class for γ < 0.2 but remains unknown for 0.20 < γ < γ c due to the rapidly increasing autocorrelation times. The authors of Ref. 58 further analyzed the dynamical behaviors of the thermal annealing method and the parallel tempering with the Metropolis update. They found that the low-temperature spin configurations obtained from thermal annealing are not always the correct state, and this might be the reason Thermodynamic results such as the specific heat provide rather limited information about the physical properties of the system at low temperatures. We further seek to understand the spin correlations or the spin structure factors that would give more important characterization of the low-temperature physical properties. These results bring us some understanding of the Ising-spin-based frustrated magnetism. In Fig. 9, we depict the (Ising) spin structure factors for different temperatures and different γ's. As long as γ > γ c = 1/4, the spin structure factors develop a reciprocal "kagomé" structure in the momentum space at low temperatures. This reciprocal "kagomé" structure can be clearly observed in Fig. 9(h) and Fig. 9(l) for γ = 0.40 and γ = 0.80, respectively. This contrasts strongly with the case of continuous spins where the reciprocal "kagomé" structure only occurs for γ = 0.5 and does not occur elsewhere (see Figs. 2, 5 and 6).
To understand the persistent reciprocal "kagomé" structure for the Ising spin with γ > γ c , we return to the geometric method in Sec. II A. The constraint in Eq. 7 cannot be satisfied for the Ising spin. To optimize the energy, one essentially minimizes S 0 + 2γS . This then requires the three spins, S 1 , S 2 , and S 3 , on the triangular corner of the unit " " to have either "two-up-onedown" or "two-down one-up" configuration, and the central spin, S 0 , to align with the minority spin. For example, for the blue triangular sublattice in Fig. 10, only the blue triangles have two-up-one-down or two-down-oneup spin configuration, not the white region. A similar requirement applies to the orange triangular sublattice. These requirements actually differ from the energy optimization for the triangular lattice Ising antiferromagnet that demands all the triangles to be either two-up-onedown or two-down-one-up. The local constraint implies that a pair of the Ising spins separated by a µ (µ = 1, 2, 3) is mostly anti-collinearly aligned. At low temperatures, the spins are fluctuating near the ground state manifold. This means that the low-temperature spin structure factor in the momentum space would be peaked at k · a µ = ±π. These peaked momenta form a reciprocal "kagomé" structure. One specific construction is given here. If one chooses k · a 1 = ±π that fixes the x components of the peaked momenta, then the Ising spins are anti-collinearly aligned along the a 1 direction. Once one a 1 -directed spin chain with anti-collinearly aligned Ising spin is realized, the two-up-one-down or two-down-oneup condition is satisfied automatically on all the triangular units boarding the spin chain. Thus, the neighboring a 1 -directed spin chains are not constrained except being anticollinearly aligned along a 1 . Thus the y component of the peaked momenta can be arbitrary. The other equivalent momenta can be generated by the crystal symmetry of the system.

IV. THERMAL EVOLUTION OF RECIPROCAL STRUCTURE FOR ISING SPINS
In the previous section, we have shown and argued that the discrete nature of the Ising moment and the frustration lead to an unexpected reciprocal "kagomé" structure in the spin structure factors at low temperatures. Moreover, this reciprocal structure is persistent with the varying of the exchange parameters. On the other hand, in several early publications, it was shown that the self-consistent Gaussian approximation seems to work even for frustrated Ising antiferromagnets [61][62][63]. This method was known to work well for frustrated XY or Heisenberg antiferromagnets. As we have found that the Ising spin behaves rather differently from the XY or Heisenberg spins at low temperatures, it is then hard to expect that the self-consistent Gaussian approximation will continue to work well for the Ising model. Apparently, if one directly applies the self-consistent Gaussian approximation to evaluate the spin structure factors, the reciprocal "kagomé" structure cannot be obtained, and the reciprocal contours would be more like the ones for the continuous spins. Then how do we reconcile the persistence of the reciprocal "kagomé" structure and the presumed applicability of the self-consistent Gaussian approximation? In addition, what is the fate of the persistent reciprocal "kagomé"s structure in the presence of the thermal fluctuations? This is addressed below.
We here fix the γ values and evaluate the spin structure factor for the Ising spins by varying the temperature. For convenience of comparison and presentation, we display the results from low temperatures to high temperatures in Fig. 11. The animations of these results can be found in the Supplemental Material [60]. We find that, at low temperatures, the spin structure factor is peaked at the reciprocal "kagomé" structure. As the temperature increases, the peak position of the spin structure factor gradually deviates from the reciprocal "kagomé" structure. In the high temperature limit, one cannot trace any signature of the reciprocal "kagomé" structure. This is expected due to the thermal fluctuations. At high temperatures, the spins fluctuate strongly, and the system deviates from the ground state manifold. Hence the stiffness of the Ising spin, which is partly responsible for the reciprocal "kagomé" structure, is conquered by thermal fluctuations.
We then perform a self-consistent Gaussian approximation to calculate the spin structure factors at finite temperatures. In this scheme, we first write down the partition function for the system, where λ i is the Lagrange multiplier that imposes the magnitude constraint for the Ising spin with |S i | = S, and S eff is the effective action and given as Here we have made a saddle point approximation by setting λ i = −β∆(T ), which is equivalent to replacing the single-spin constraint with a global spin constraint. The uniform saddle point λ is based on the property of the homogeneous paramagnetic state that respects all the lattice symmetry. Thus the spin correlation function is then given as where I 2×2 is a 2 × 2 identity matrix. The saddle point parameter ∆(T ) is determined self-consistently from the saddle point equation, Clearly, in this scheme, the momentum information arises from the exchange interaction matrix J (k). Thus, the spin structure factor is mostly weighted around the degenerate contours of the lowest eigenvalues of J (k). This is a bit analogous to the spiral spin liquid regime for the continuous spins in Sec. II. In fact, we have computed the spin structure factors within the self-consistent Gaussian approximation for the continuous spins in Fig. 5. Our results from the Ising spins are depicted in Fig. 11. In the upper (lower) set of panels of Fig. 11, we choose γ = 0.40 (γ = 6.0) such that the contour is around the center (corner) of the Brillouin zone. The comparison with the numerical results is better at the high temperatures and is poorer at the low temperatures. This behavior is expected. In the high temperatures, the Ising spins are strongly fluctuating thermally, and the spin magnitude constraint does not play a strong role. The system widely leaves the ground-state manifold such that the reciprocal "kagomé" structure is no longer visible, but the correlations between the spins are still preserved from the exchange interaction part. The self-consistent Gaussian approximation captures this high temperature correlated regime and gives a qualitatively reasonable description of the spin structure factors. As a comparison for the continuous spins in Fig. 5, the spiral contours are well captured by the self-consistent Gaussian approximation from the high temperatures to the low temperatures except that the intensity is not well obtained in the lowtemperature limit.

V. DISCUSSION
We have shown that the spiral spin liquid regime holds for both XY and Heisenberg spins in the frustrated regime at low but finite temperatures. We further identified the finite-temperature chiral transition for both cases. For the Ising spin, a reciprocal "kagomé" structure emerges in the low-temperature spin structure factors in the momentum space, and persists for a range of exchanges in the frustrated regime which resembles the spiral spin liquid regime. This is understood from the stiffness of the discrete Ising spins and the frustration due to the competing interactions on the honeycomb lattice. Moreover, the reciprocal structure evolves from the "kagomé" structure into the contours demanded by the soft spin analysis as the temperature is increased.
It is readily evident as well as illuminating to notice that, if one adds an anisotropic term such as D(S x i ) 2 [or D(S z i ) 2 ] to the model for the XY (or Heisenberg) spin, the system will behave like the Ising spin with an increasing and positive D. Thus one anticipates that the spiral contour of the XY and Heisenberg cases will crossover to the reciprocal kagomé structure. The details of such a crossover process were not pursued in this paper, and we expect that it could be observed in real materials of the spiral spin liquid regime through the tuning of such spin anisotropy, e.g. via hydrostatic pressure [9][10][11].

A. Experimental connection
We here discuss the experimental connection. As far as the materials' relevance, many honeycomb lattice van der Waals magnets have been proposed and synthesized. For a recent review of them, one can refer to Ref. [1]. More recently, rare-earth honeycomb lattice magnets were proposed and studied [45][46][47]. For the rare-earth magnets, the exchange interaction is often anisotropic and short ranged. The exchange coupling decays rapidly with the distance and often becomes quite weak beyond the first neighbor. Thus, the frustration on the rare-earth honeycomb lattice magnets often arises from the anisotropic exchange, rather than the competing further-neighbor interactions. From this perspective, it seems a bit difficult to expect the physics in this paper to be directly relevant for these honeycomb lattice rare-earth magnets. In fact, the Yb honeycomb lattice magnet YbCl 3 was found to have a simple Néel state [64], indicating the dominance of the first-neighbor interaction [65]. However this is not the end of the story for the honeycomb lattice rare-earth magnets.
A series of vdW rare-earth chalcohalides were recently synthesized [48]. These materials are not of planar honeycomb structure like YbCl 3 [64]; instead, they are formed by the bilayer triangular lattice of the rare-earth moments. The bilayer triangular lattice is A-B stacked and is equivalent to a honeycomb lattice. The firstneighbor and second-neighbor distances are close. For HoOF, the first neighbor is 3.57Å, and the second neighbor is 3.80Å. Some of these compounds certainly support highly anisotropic spin interactions such as the Kitaev interaction [45][46][47]; the Ho-based, Dy-based, and even Er-based ones could actually provide the Ising local moments [48] as the effective moments of these ones are close to the fully polarized atomic values. If we assume that both the first-and second-neighbor Ising interactions are due to the dipole-dipole interaction, this will put these materials in the frustrated regime with γ > γ c . Thus, if the remaining weak interactions beyond the first and second neighbor Ising interactions do not play a significant role, we expect our prediction including the reciprocal "kagomé" structure and the thermal evolution of the reciprocal structure in this paper to be applied to these compounds. Experimentally, this requires an inelastic neutron scattering measurement of the dynamic spin correlation for a range of energies that is integrated to yield the equal-time spin correlation.
In this family of materials, the Yb and Sm ones seem to be more quantum, and most likely to carry quantum spin-1/2 Kramers doublets. According to a more microscopic study of the exchange paths [49], the exchange interaction between the Yb local moments may sometimes be Heisenberg-like. In that case, the Yb compound may realize a competing spin-1/2 J 1 -J 2 Heisenberg model, and the Dzyaloshinskii-Moriya interaction could play some role here. The spiral spin liquid regime may apply to the finite temperature regime of the spin-1/2 J 1 -J 2 Heisen-berg model, but the ground state of this model can be interesting as well [66]. In addition, in the J 1 -J 2 model with the Heisenberg spin, adding an external field perpendicular to the lattice could bring about a skyrmion crystal [67,68], and thermal Hall effect of magnons could be generated. The magnetic skyrmion was suggested to have a great stability on account of its topological protection, which makes it an excellent information carrier, and its topological properties induce a variety of emergent behaviors attracting research interest. In the vdW materials, this will not only have theoretical value but also promote the development of device applications.
The Co-based honeycomb lattice antiferromagnets were recently proposed to be candidates for Kitaev materials [69,70] due to the spin-orbit-entangled local moment. In fact, the Co-ions were also known to support Ising interactions. This has been found in the quasione-dimensional Ising magnets CoNb 2 O 6 , BaCo 2 V 2 O 8 and CaCo 2 V 2 O 8 [71][72][73][74]. It would also be interesting to search for Co-based honeycomb lattice antiferromagnets that realize Ising spins.
Another set of compounds that are not really magnets and do not rely on the spin degrees of freedom are the Fe-based mixed valence compounds. The system contains both Fe 2+ ions and Fe 3+ ions, which can be equivalently treated as Ising spin, and thus realizes an interacting Ising spin system. This idea has been applied to understand the electron charge physics in LuFe 2 O 4 and YbFe 2 O 4 [61,62]. The lattice in these two compounds is not honeycomb but frustrated. If an Fe-based mixed valence compound with a honeycomb lattice is found, our result can be applied, and the spin correlation should be replaced by the electron density-density correlation, which can be detected by X-ray scattering.

B. Frustrated Ising magnets
Here, we discuss the implications of our results for frustrated Ising magnets. In this paper, we show that the high-temperature spiral spin liquid regime is not quite sensitive to the spin dimensions. For both the continuous spin and the Ising spin, the same kind of momentumspace contours demanded by the exchange interaction appear in the spin structure factor. This is well captured by the self-consistent Gaussian approximation. However, at low temperatures, the Ising spin in the honeycomb lattice of our problem develops an unconventional reciprocal "kagomé" structure. This means that frustrated Ising magnets may contain more interacting correlated spin structures at low temperatures than the ones expected from the self-consistent Gaussian approximation. Frustrated Ising magnets in an A-B stacking multilayer triangular lattice or an A-B-C stacking multilayer triangular lattice, which were studied in LuFe 2 O 4 and YbFe 2 O 4 [61,62], may contain some new ingredients at low temperatures beyond what has been found using the high-temperature Gaussian approximation, and may be worth a careful examination. A similar situation may also occur for the diamond lattice Ising antiferromagnet in the frustrated regime.