Phase Transitions for Cuboc Orders in Stacked Kagome Heisenberg Systems

Using the event-chain Monte Carlo (MC) algorithm, we investigate phase transitions of the stacked Kagome Heisenberg systems with classical vector spins including up to the 3rd-nearest-neighbor couplings. In particular, we focus on two types of non-coplanar spin orders ---cuboc1 and cuboc2 orders, both of which have twelve-sublattice structures accompanying the translational-symmetry breaking. We perform event-chain MC simulations up to L=72, where L represents the linear dimension of the stacked Kagome lattice, and then find that the cuboc1 transition shows 2nd-order-transition behaviors with a tendency to a weak-1st-order transition up to L=72, while the cuboc2 transition is basically described by the 1st-order transition. We then discuss the above transitions in connection with the effective Landau-Ginzburg-Wilson theory with the O(3)xO(3) symmetry.

Frustrated spin systems often provide fascinating physics associated with nontrivial spin orders having non-colinear or non-coplanar spin structures. Recently, a variety of exotic ground-state spin orders were classified for the Kagomelattice classical Heisenberg antiferromagnets including up to the 3rd-nearest-neighbor couplings. 1) Among various spin orders for the Kagome Heisenberg systems, we particularly focus on the cuboc orders, which are defined as non-coplanar spin configurations with twelve sublattice structure specified by the triple-q structures in the momentum space. 2,3) Moreover, experiments on CsCrF 4 , 4, 5) NaBa 2 Mn 3 F 11 6) and Cu 3 Zn(OH) 6 Cl 2 7, 8) suggest the relevance of the cuboc orders to realistic experimental situations.
Although the planar-lattice spin models with continuous symmetries have no finite-temperature order, the inter-layer coupling may lift the exotic ground-state orders into the finitetemperature orders with the spontaneous symmetry breaking, which provide rich physics associated with the phase transitions and critical behaviors. Actually, the stacked triangularlattice Heisenberg antiferromagnets have been extensively studied so far, where the effective Landau-Ginzburg-Wilson (LGW) theory based on the O(3)×O(2) symmetry plays a central role to understand the universality reflecting the double-q structure of the planar 120 • order. [9][10][11][12][13][14][15][16][17][18][19][20][21] However, phase transitions associated with the non-coplanar spin orders, which are characterized by the triple-q structures, 11) have been less explored, because systematical analysis of the complicated non-coplanar orders is difficult from the numerical simulation point of view.
A crucial point on the stacked Kagome systems is that they can realize two similar but different types of cuboc orders, which are mentioned as cuboc1 7) and cuboc2 2) in Ref. 1; although these two cuboc orders have the same symmetry in the spin space, the triple-q wave vectors specifying their Bragg peak points are located at different positions in the momentum space. Thus, the stacked Kagome systems can be an interesting stage for investigating the universality of the cuboc transitions, which may be described by the O(3)×O (3) LGW theory. In addition, to clarify the stability of the orders would be essential for experiments on the Kagome lattice compounds.
In this letter, we investigate phase transitions of the stacked  Kagome Heisenberg systems depicted in Fig. 1(a), using the event-chain Monte Carlo (MC) algorithm, which was originally introduced for hard core particle systems 22,23) and recently developed to the classical spin systems. 24,25) The eventchain method is a rejection free MC algorithm satisfying the global-balance condition, and enable us to achieve the proper relaxation even for the complicated spin orders in frustrated spin systems. We perform MC simulations for the stacked Kagome systems exhibiting the cuboc orders up to about 10 6 spins, and then analyze the natures of the cuboc transitions in detail. We also discuss the connection to the O(3)×O (3) LGW theory. The planar structure in the stacked Kagome systems is depicted in Fig. 1(a), where J 1 is the nearest-neighbor coupling and J 2 denotes the half of the next-nearest-neighbor couplings. We also introduce the 3rd-nearest-neighbor couplings J d that run along the diagonal directions of the hexagons in the Kagome lattice. Then, the Hamiltonian is explicitly written as (Color online) (a) Cuboc1 spin structure: the four tilting 120 • planes at the small numbered triangles in the Kagome lattice form the tetrahedron in Fig. 1(b), accompanying the translational symmetry breaking. In the extended Brillouin zone corresponding to the outer hexagon in the upper right panel, the cuboc1 is represented as the modes indicated by the solid symbols. (b) Cuboc2 spin structure, the 120 • planes are located on the triangles of the J 2 couplings. In the Brillouin zone of the inner hexagon, the cuboc2 is represented as the solid symbols at the M points. The numbers of the triangles correspond to those of the surface triangles in Fig. 1 where S denotes the vector spin of the O(3) symmetry with |S| = 1 and J c represents the inter-layer couplings. We basically assume |J c | = 1. Note that J c causes no frustration effect, implying that spin structures of the ordered phases are determined by the frustrating couplings in the Kagome plane.
In order to see the formation of the cuboc orders for Eq. (1), we consider the J 1 − J d Kagome system with J 2 = 0 for a while. As shown in Ref. 1, we can see that the cuboc1 is observed for J 1 > 0 and J d > 0, where the J d couplings select the cuboc1 state from the highly-degenerating ground states of the Kagome system having only the nearest-neighboring J 1 couplings. Figure 2(a) illustrates the spin configurations of the cuboc1 order in the real space, where the J d couplings stabilize the staggered spin arrays along the diagonal directions. In order to reduce the energy due to the frustrating J 1 couplings on the small colored triangles in Fig. 2(a), the spin orientations in the three staggered arrays relatively tilt to form the 120 • structure, accompanying the translational symmetry breaking. Then, the four tilting 120 • planes on the small triangles labeled by 1 to 4 in Fig. 2(a) can be mapped into the tetrahedron in Fig. 1(b). As shown in the extended Brillouin zone in the upper right panel of Fig. 2(a), the cuboc1 configuration is characterized by the triple q located at the internally dividing point in the ratio 3:1 between the Γ and K points. 1) In the region −2J d < J 1 < 0 with J d > 0, on the other hand, we find that the cuboc2 is realized as the ground state. Although J 1 is ferromagnetic in this case, the three staggered spin arrays along J d directions are frustrating on the hexagons and the spins on the large colored triangles in Fig. 2(b) form the 120 • structure. Then, the four 120 • planes labeled by 1 to 4 in Fig. 2(b) form the cuboc2 order accompanying translational symmetry breaking, which is represented as the same tetrahedron in Fig. 1(b). The difference between the cuboc1 and the cuboc2 is in the real-space distribution of the 120 • planes. In contrast to the cuboc1 case, the cuboc2 is described by the triple q corresponding to the M points in the Brillouin zone of the Kagome lattice, as shown in the upper right panel of Fig. 2(b). In this sense, it is an interesting problem to check the universality of the two types of cuboc transitions.
Another interesting point on the above cuboc2 phase is that, if the J 2 coupling turns on and J d decreases, it can be adiabatically connected to that in the Kagome-triangular spin system, where the positive J 2 coupling strongly stabilizes the 120 • structure on the J 2 triangles. 5,26) Moreover, the phase transitions in the Kagome-triangular system were intensively investigated in the context of the coupled spin tubes by MC simulations up to about 10 5 spins, 5) where the crossover depending on the ratio |J 1 |/J 2 between the weak-1st-order and 2nd-order transitions was reported. However, the previous MC results still contained a significant finite-size effect, so that a large scale MC simulation was required for verifying the universality of the cuboc2 transition.
In order to discuss the phase transitions for Eq. (1), we set up three order parameters characterizing the cuboc orders. The sublattice magnetization M is defined as the expectation value of spins respectively for the 12 sublattices. The sublattice vector spin chirality κ V is defined for each triangle having the 120 • spin structure as κ V ≡ where A, B, and C are the index of the spins in the 120 • -structure triangle. As seen in Fig. 2, the four triangles are contained in the magnetic unit cell having 12 sublattices. Also, the sublattice scalar spin chirality κ S , which detects the non-coplanar spin structure, is defined as κ S ≡ S A · [S B × S C ], where S A , S B and S C should belong to the three adjacent surface triangles of the tetrahedron in Fig. 1(b). For the cuboc1, κ S is defined on the large triangles embedded in the hexagons in Fig. 2(a), while for the cuboc2, it is located on the small triangles in Fig. 2(b). The alternating locations of the triangles of κ V and κ S can be viewed as a duality relation between the cuboc1 and cuboc2.
In order to analyze the finite temperature transitions for the cuboc orders in the stacked Kagome Heisenberg systems, we employ the event-chain algorithm [22][23][24][25] combined with the parallel tempering, 27) which enable us to achieve the relaxation to the thermal equilibrium even for the highly frustrated systems. We then calculate the susceptibilities of the order parameters, χ M = |M| 2 /T , and χ V = |κ V | 2 /T , where T is a temperature. 28) In actual simulations, the typical number of the sampling is about 10 5 after the initial relaxation of 10 4 MC steps, and the maximum system size is L = 72, where L is the linear dimension of the stacked Kagome system in the unit of the spin triangle (unit cell of the Kagome lattice). Note that we terminate a single event-chain update for a spin rotation axis randomly chosen from the x, y, z directions, if the total rotation angle reaches 6πL 3 /100.
Let us first discuss the cuboc1 transition for the J 1 − J d stacked Kagome system with J 2 = 0 and J c = −1 fixed. In  Fig. 3(a), which basically indicates the 2nd-order transition. Indeed, the transition temperature and the critical exponents are extracted as T c = 2.1549(3), ν = 0.63(2), and γ = 1.22 (5), with use of the Bayesian approach. 29) Taking account of the L-dependence of the estimated exponents, we adopt ν = 0.60(4) and γ = 1.17 (10) (2) as the exponents for J 1 = 1, J d = 3.2 and J c = −1. In Fig.  3(b), we also present the FSS plot for χ V , which yields the exponents ν V = 0.59(2) and γ V = 0.56(3) at the same T c within the numerical accuracy. 28) Assuming the 2nd-order transition, moreover, we have verified that Eq. (2) consistently reproduces the FSS plots in 3.2 ≤ J d ≤ 100. The results are summarized in the phase diagram of Fig. 4(a). Note that, in the region 0 < J d 1.0, the MC simulation fails in the relaxation to the thermal equilibrium state, due to the significant frustration originating from the J 1 Kagome couplings.
Although the above FSS scaling analysis basically suggests the 2nd-order transition for the cuboc1, we should remark the possibility of a weak-1st-order transition. This is because, for Eq. (2), the anomalous dimension η ≡ 2 − γ/ν can be negative within the error bars, which could be a signature of the 1storder transition. Also, the slightly large error bars in Eq. (2) originate from a non-negligible weak J d dependence involved in the FSS analysis within L = 72, which might suggest a crossover to the 1st-order transition. In order to establish the cuboc1 transition, we need a larger-scale simulation, which is an important future issue.
We next turn to the cuboc2 transition for the J 1 − J d stacked Kagome systems in the ferromagnetic J 1 regime. In 0 ≤ −J d /J 1 ≤ 0.5, the ground state is ferromagnetically ordered, while the cuboc2 order is realized for −J d /J 1 > 0.5. Thus, the finite-temperature transition to the cuboc2 order can be observed in −J d /J 1 0.5. In the following, we fix J 1 = J c = −1 and vary J d . The phase diagram determined by the MC simulations up to L = 72 is summarized in Fig. 4(b). In J d 0.5, we have observed the ferromagnetic transition, which is consistent with the ground-state behavior. Note that this ferromagnetic transition is verified to be in the 3D ferro- magnetic Heisenberg universality. For J d 0.5, on the other hand, we find that the cuboc2 transition actually appears at a finite temperature. Let us discuss the nature of the cuboc2 transition in details. An important point is that the 1st-order transition is confirmed in the small J d regime (0.5 J d ≤ 0.8). Figure 5(a) illustrates the energy histogram for J d = 0.7 at T = 0.9147, where the double-peak structure emerges for L = 72. In Fig.5(b), we can also confirm the double peaks of the scalar-spin-chrality histogram at κ S ≃ 0 and 0.013. 30) Taking account of the temperature range where the double peaks are observed, we adopt T c = 0.9147(2) as a transition temperature. As J d increases beyond J d ∼ 0.8, the transition crossovers to the 2nd-order transition like behaviors. However, the FSS plots of χ M for J d > 0.8 yield non-universal exponents for χ M gradually gliding in ν ≃ 0.51 − 0.72 and γ ≃ 0.88 − 1.38. Moreover, we find that the FSS plot of χ V fails in the range J d > 0.8. On the basis of these nonuniversal behaviors, we think that the cuboc2 transition is basically described by the 1st-order transition in the bulk limit, though larger-scale simulations are required for the direct evidence of the 1st-order transition.
As mentioned before, the cuboc2 order in the J 1 − J d stacked Kagome system is adiabatically connected to that in the Kagome-triangular system with J 1 < 0 and J 2 > 0. In the following, moreover, we assume the antiferromagnetic interlayer coupling J c = 1, for the purpose of a direct comparison with the coupled-spin-tube system previously studied in Ref.5, where it was reported that the 2nd-order-transition behaviors in −J 1 /J 2 0.8 crossover to the 1st-order transition in 0.8 −J 1 /J 2 1.0 within L ≤ 36. Note that the interlayer coupling cause no frustration and thus, the nature of the transitions is irrelevant to the sign of J c .
As L increases up to L = 72, in this paper, we have found that the FSS plot fails in extracting the critical exponents even for −J 1 /J 2 = 0.5 and then the energy histogram exhibits the double-peak structure at T c = 0.6470(1) estimated, as shown in Fig. 6(b). Also we have confirmed that the histograms of the spin chiralities exhibit the double-peak structures at the same T c . These results imply that the region where the 1storder transition appears is extended to the small −J 1 /J 2 side, as L increases. We then illustrate the T − J 1 phase diagram for the stacked Kagome-triangular system of J 2 = 1 and J c = 1 in Fig. 6(a), where the 1st-order transition is confirmed in 0.3 ≤ −J 1 /J 2 ≤ 1.0. Note that the double-peak structures of the histograms can be also confirmed in a wide range of the inter-layer coupling J c within the relatively small system sizes compared with the J 1 − J d Kagome case. A reason for this weak finite-size effect is that the J 2 coupling strongly stabilizes the 120 • structure, so that the effective length scale at T c could be shorter than that of the J 1 − J d Kagome systems.
In this letter, we have investigated the phase transitions in the stacked J 1 − J d Kagome Heisenberg systems, which exhibit the cuboc1 order for J 1 > 0 and the cuboc2 order for −2J d J 1 < 0. We also analyzed the cuboc2 transition in the stacked Kagome-triangular spin systems. An important point is that the cuboc1 and cuboc2 orders are represented as the triple-q structures in the momentum space, though they have similar but distinct real-space configurations of the noncoplanar spin structure. From the effective field theoretical view point, then, the phase transitions are basically described by the O(3) × O (3) LGW theory, 11,14) for which the 1st-order transition is expected on the basis of the 4 − ε expansion. 20) In this sense, the 1st-order transition for the cuboc2 order, which is supported by the double-peak structure of the energy histogram at the transition temperature, is consistent with the field theoretical analysis. Although the 2nd-order-like behaviors are observed for the large J d region of the J 1 − J d Kagome system within L = 72, we think that the nonuniversal exponents extracted is a signature of the 1st-order transition.
In contrast, the FSS analysis for the cuboc1 transition suggests the 2nd-order transition behaviors within L = 72, which might contradict to the analysis of the O(3) × O (3) LGW theory. However, the resulting anomalous dimension η = 2 − γ/ν estimated for the cuboc1 transition can be negative within the error bars, which could also be a signature of the 1st-order transition. We think that the discrepancy of the finite-size effects between the cuboc1 and cuboc2 transitions are basically attributed to the nonuniversal features depending on the model parameters. For the Kagome-triangular cases where J 2 coupling strongly stabilizes the 120 • structure embedded in the cuboc2 order, we have actually observed the clear evidence of the 1st-order transition. In order to confirm the nature of the cuboc1 transition, larger scale numerical simulations are clearly required. Also, it is an important problem to investigate the universality/nonuniversality for other non-coplanar spin orders, 31) taking into account the various field-theoretical approaches. 12,15,16,19) Although the phase transitions and the critical phenomena are a long standing issues in physics, the quantitative analysis for the non-coplanar spin orders has been less explored, since the targeting systems are basically highly frustrated. We have demonstrated that the event-chain MC simulation is a very powerful tool for the quantitative analysis of the stacked Kagome Heisenberg systems, which are typical systems exhibiting the phase transitions to nontrivial non-coplanar spin orders. Our results stimulate further researches of the phase transitions associated non-coplanar spin orders, which involve a new frontier of statistical mechanics and condensed matter physics.