Finite-temperature properties of the Kitaev-Heisenberg models on kagome and triangular lattices studied by improved finite-temperature Lanczos methods

Frustrated quantum spin systems such as the Heisenberg and Kitaev models on various lattices, have been known to exhibit various exotic properties not only at zero temperature but also for finite temperatures. Inspired by the remarkable development of the quantum frustrated spin systems in recent years, we investigate the finite-temperature properties of the S = 1/2 Kitaev-Heisenberg models on kagome and triangular lattices by means of finite-temperature Lanczos methods with improved accuracy. In both lattices, multiple peaks are confirmed in the specific heat. To find the origin of the multiple peaks, we calculate the static spin structure factor. The origin of the high-temperature peak of the specific heat is attributed to a crossover from the paramagnetic state to a short-range ordered state whose static spin structure factor has zigzagor linear-intensity distributions in the momentum space. In the triangular Kitaev model, the “order-by-disorder” due to quantum fluctuation occurs. On the other hand, in the kagome Kitaev model it does not occur even with both quantum and thermal fluctuations.

The S = 1/2 Kitaev model on the honeycomb lattice (HL) dose not have geometric frustration but has frustration effect arising from the bond-dependent Ising interactions 54 , called exchange frustration. In this model, the S = 1/2 spins are divided into localized Majorana fermions composing Z 2 fluxes and itinerant Majorana fermions [55][56][57] . Its ground state exhibits an exact quantum spin liquid with topological order. At finite temperatures, there is a distinct double peak in the specific heat 58 . The origin of this double peak is described below: the high-temperature peak is caused by freezing the itinerant Majorana fermions and the low-temperature peak is caused by freezing the localized Majorana fermions 58 . Because of clear difference of their energy scales, a 1/2plateau-like anomaly appears in the temperature dependence of the entropy. This phenomenon corresponds to a fractional excitation of the spins. Moreover, such a phenomenon has been found even in the spin S > 1/2 and mixed spin systems 59,60 , even though the spin degree of freedom cannot be decomposed into Majorana fermions. Furthermore, finite-temperature properties of the Kitaev-Heisenberg (KH) model have also been studied on the HL 61,62 .
The S = 1/2 KH models on the KL and TL, having both the geometric frustration and exchange frustration, have been studied mainly for the ground state [63][64][65][66][67][68][69] . In the KL-KH system, it has been proposed that there are two quantum spin-liquids, a canted ferromagnetic, and a q = 0, 120 • ordered phases 63 . Whereas in the triangular KH system, it has been proposed that there are Z 2 vortex crystal, nematic, dual-Z 2 vortex crystal, ferromagnetic, and dual-ferromagnetic phases [65][66][67][68][69] . However, finite-temperature properties in the KH models on the KL and TL have hardly been investigated. There is possibility that multiple peaks in the temperature dependence of the specific heat and new crossover phenomena exist, because such phenomena have been confirmed in the HL-Kitaev and KL-Heisenberg models. Therefore, it is important to investigate the finite-temperature properties of these models.
The finite-temperature Lanczos method (FTLM) is a useful technique for calculating finite-temperature properties 70,71 . However, this method has a problem that the accuracy becomes worse at low temperatures 71 . Therefore, we need to overcome this problem. In this paper, we first propose two methods to improve the FTLM. We name the methods replaced finitetemperature Lanczos method (RFTLM) and orthogonalized finite-temperature Lanczos method (OFTLM). Using these improved FTLMs, we next calculate the specific heat, entropy, and static spin structure factor (SSSF) to investigate the finite-temperature properties of the S = 1/2 KH model on the KL and TL. In the kagome system, the specific heat exhibits multiple-peak structures at finite temperatures for 0 ≤ θ ≤ 0.5π, where θ = arctan(K/J) with J (K) being the Heisenberg (Kitaev) interaction. To clarify the origin of the multiplepeak structure of the specific heat, we analyze the SSSF at finite temperatures for the N = 36 cluster using the RFTLM. From the analyses, we find that the highesttemperature peak of the specific heat for 0.1π ≤ θ ≤ 0.4π is originated by a crossover from the paramagnetic state to a state whose SSSF intensity shows direction distribu- tion in the momentum space. On the other hand, one of low-temperature peaks for 0.1π ≤ θ ≤ 0.4π is expected to be a signature of the emergence of a q = 0, 120 • order state. However, at θ = 0.5π (Kitaev limit), the q = 0, 120 • order does not appear. In the triangular system, we find that there is a double-peak structure in the specific heat for 0.25π ≤ θ ≤ 0.5π. The origin of the doublepeak structure is the same as the kagome system. At θ = 0.5π, the ground state exhibits a stripe order due to the "order-by-disorder" mechanism unlike the kagome system.
The arrangement of this paper is as follows. In Sec. II, we describe our S = 1/2 KH models on KL and TL. In Sec. III, we first explain the standard FTLM, then we explain the RFTLM and OFTLM developed by us. In Sec. IV, the results of the specific heat, entropy, and SSSF for the KL and TL are shown. In Sec. V, we discuss difference between the honeycomb, kagome, and triangular systems for the origin of the multiple-peak structures in the specific heat and we focus on characteristic of the Kitaev model on the KL. Finally, a summary is given in Sec. VI.

II. MODEL
The Hamiltonian of the KH model is given by where S i is a quantum spin operator with S = 1/2 at site i. J i,j represents the nearest neighbor interactions as shown in Fig. 1(up panel) for the KL and Fig. 1(down panel) for the TL. J i,j takes one of the three anisotropic interactions, J X = diag(J + K, J, J) (yellow bonds), J Y = diag(J, J + K, J) (light green bonds), and J Z = diag(J, J, J + K) (blue bonds), where K and J correspond to the energy of the Kitaev and Heisenberg interactions, respectively. We introduce the parametrization (J, K) = (I cos θ, I sin θ), where I is the energy unit (I = 1). In the present study, we focus on 0 ≤ θ ≤ 0.5π.

A. Finite-temperature Lanczos method
In this section, we describe the standard FTLM 70,71 . The FTLM has been used to study the finite-temperature properties of various lattice models 48,51,[72][73][74][75][76][77][78] . The partition function Z(T ) of the canonical ensemble at temperature T is expressed as follows: where N st is the dimension of H, |n is an arbitrary normalized vector, β is the inverse temperature 1/T (k B = 1), E i is an eigenenergy of H, |Ψ ik is an eigenvector with E i , d i is a degree of degeneracy of the state with E i , and N d represents the number of the eigenenergies, which satisfies N st = The FTLM introduces two approximations for (2). The first one is to replace the summation of n with random sampling r with R times. The second one is for the summations of i and k. Both the summations are replaced by the Krylov subspace with dimension M . In the FTLM, the partition function and general operator A are approximated as follows: where |V r is a normalised random initial vector and |ψ r j ( (r) j ) are an eigenvector (eigenvalue) in the M -th Krylov subspace for H. We note that |V r is formally given by |V r = N d −1 i=0 di k=1 η rik |Ψ ik using the exact eigenstate |Ψ ik , where η rik is a random value that satisfies di k=1 |η rik | 2 = 1 for the normalization. For the energy E(T ), specific heat C(T ), and entropy S(T ), the following general expressions are useful: . From these equations, E(T ) and C(T ) calculated by the FTLM are given by At high temperatures, a few sampling R is enough for obtaining high accuracy since the error of all physical quantities is Equation (8) does not give an exact value, and if A is noncommutative with Hamiltonian such as the SSSF, Eq. (9) also does not give an exact value. These errors are expected to be O(1/ √ R) 71 . Therefore, a very large number of samplings are required to obtain good accuracy at low temperatures. The low-temperature Lanczos method 79 is known as one of the solutions to this problem. However, this method has a difficulty for large-scale calculations because it requires huge random access memory to keep all vectors in the Krylov subspace with M . Therefore, we try to improve the accuracy of the FTLM at low temperature in two ways: the RFTLM and OFTLM.

B. Replaced finite-temperature Lanczos method (RFTLM)
In the standard Lanczos method, we can obtain several low-lying eigenstates with N E levels whose energy is given by (r) i (i = 0, 1, · · · , N E − 1), but cannot judge degeneracy of each level. Therefore, N E −1 and each eigenvector should be written generally Here, we assume that the obtained energy (r) i is independent of sampling r, i.e, (r) i = E i , although the corresponding eigenvector may depend on the sampling |ψ r i = |Ψ r i due to possible degeneracy. Then we can rewrite the expression (4) as follows: Comparing the first term on the right hand side of Eq. (10) with Eq. (3), we come up with replacing The replacement (11) leads to the partition function of the RFTLM The first term in Eq. (12) is the same as the exact partition function Z(T ) (3), for i < N E . This indicates that Z(T ) RFTL (12) is more accurate than Z(T ) FTL (4). In a similar way, A (T ) FTL can be improved in accuracy by replacing We can obtain the exact eigenstates |Ψ ik with E i by the several kinds of exact diagonalization (ED) methods such as the thick-restart Lanczos method 80 , band-Lanczos method 81 , locally optimal block preconditioned conjugate gradient method 82 , and root-shifting method 83 .
By performing the RFTLM, S(T → 0) RFTL and A (T → 0) RFTL become an exact value ln(d 0 ) and d0 k=1 Ψ 0k |A|Ψ 0k /d 0 , respectively. Therefore, accuracy at low temperatures using the RFTLM would be extremely improved as compared with the standard FTLM. The efficacy of the RFTLM is confirmed in Sec. III D. However, in the RFTLM, it is necessary to know the degeneracy d i in order to perform the summation of i in Eqs. (12) and (14). We also should be careful about pseudo-eigenvalues so-called "ghost" eigenvalues caused by the presence of the machine epsilon. If there are the ghost eigenvalues, it is necessary to change N E in the second term of Eqs. (12) and (14) to N E + N g , where N g is the number of the ghost eigenvalues less than E N E −1 . We develop a new method in the next section to overcome these problems.

C. Orthogonalized finite-temperature Lanczos method (OFTLM)
In this subsection, for simplicity, we include the index k for degeneracy into the index i hereafter; rewriting η rik ⇒ η ri and |Ψ ik ⇒ |Ψ i . Thus the random vector |V r reads |V r = Nst−1 i=0 η ri |Ψ i . In the OFTLM, we first calculate several low-lying exact eigenvectors |Ψ i with N V levels (E 0 ≤ E 1 ≤ · · · ≤ E N V −1 ) before performing the FTLM. We next use following modulated random vector: with normalization Here, |V r is orthogonal to the states |Ψ i for i < N V . Therefore, FTLM using |V r as the initial vector is equivalent to applying the method to a Hilbert space excluding |Ψ i through and A (T ) of the OFTLM are obtained by adding exact values coming from |Ψ i to the FTLM result obtained by using |V r as an initial vector Since |Ψ i obtained by the ED methods would be slightly different from the exact vectors because of the machine epsilon, some of (r) j in FTLM using |V r may become, for example, E 0 , which should not appear. In practical use, this is no problem since | V r |ψ r j | for such a E 0 becomes extremely small (∼ machine epsilon). We can see that Z(T ) OFTL and A (T ) OFTL are close to the exact values at low temperatures. We emphasize that in the OFTLM we do not need to know degeneracy d i in |Ψ i and can make M smaller compared to the FTLM and RFTLM. The efficacy of the OFTLM is confirmed in Sec. III D.

D. Confirming the efficacy of the RFTLM and OFTLM
We perform benchmark calculations for the standard FTLM, RFTLM, and OFTLM. We calculate S(T ), C(T ), and the z component of SSSF, with the position vector r j and r k . All FTLMs are performed with M = 90 and R = 50. Here, we note that M = 90 is large enough to obtain the ground state of the N = 12 kagome system. Calculated results are shown in Fig. 2. The standard errors of the FTLMs using the jackknife technique 84 are represented by blue shaded regions in Fig. 2. We can see that the accuracies of the RFTLM and OFTLM are clearly better than that of the standard FTLM for all physical quantities. Therefore, we succeed in improving the FTLM. Furthermore, we compare the standard FTLM and OFTLM in detail using S z q (T ) in Fig. 3. In the standard FTLM, the accuracy for M = 30 is very poor at low temperatures as shown in the left panel of Fig. 3 because of small M that is not enough to make a convergence to the ground state. On the other hand, high-precision results can be achieved in the OFTLM even for the same M (see the middle panel of Fig. 3), since the contribution from low-energy sectors are added separately as shown in Eqs. (17) and (18). For this reason, the OFTLM give a good convergence quicker then other FTLMs. In the OFTLM with larger M , the eigenvalues less than E N V −1 and the ghost eigenvalues appearing in the first terms of Eqs. (17) and (18) may affect S z q (T ). In order to investigate these effects, we also perform the OFTLM with very large M = 5000 (> N st ). We can see that there is no effect on S z q (T ) as shown in the right panel of Fig. 3. This means that the OFTLM is not only a highly accurate method but also a user-friendly method because one can choose M without checking the convergence of eigenvalues in each Lanczos sampling.

A. Conditions of numerical calculation
In present study, we calculate C(T ), S(T ), and S z q (T ) using the RFTLM for N = 36 and the OFTLM for N = 24 and N = 30. The N = 24, N = 30, and N = 36 clusters are shown in Fig. 1 for the KL and TL. Finite size effects can be reduced by using large size and highly symmetric clusters such as N = 36. We emphasize that the improved FTLMs with high accuracy make finite-size effects at low temperatures very clear. To calculate the excited states required for using the improved FTLMs, we use the restarted Lanczos method with the root-shifting method. Table I shows detailed conditions for improved FTLM calculations.
For large clusters such as N = 36, it is time-consuming to prepare several eigenvectors with N E > 1 or N V > 1. Furthermore, one has to be careful the appearance of the ghost eigenvalues in such a huge calculation. To avoid these difficulties, we decide to use the RFTLM with N E = 1, where we need to calculate the ground state only. The accuracy of the N E = 1 result will be confirmed in the next section.

B. Kagome lattice
We first discuss the efficiency of the RFTLM for the N = 36 kagome system at θ = 0.2π. Figure 4 shows S z q (T ) at q = (2π, 2π/ √ 3) using the standard FTLM and RFTLM. In the standard FTLM, there is large error at low temperatures and an average value at T = 0 deviates from the exact one. On the other hand, in the RTLM, error bars become less than the width of the line for all temperatures and an average value converges to the exact one at T = 0. This clearly demonstrate that our improved FTLMs work well even for the N = 36 system. We emphasize that the error of the FTLMs becomes almost less than the line width in all the results shown below. Figure 5 shows the calculated results of C(T ) (left panels) and S(T ) (right panels) for 0 ≤ θ ≤ 0.5π at N = 24, 30, and 36. C(T ) exhibits the multiple-peak structures in all θ and N . For T > 0.2 and all θ, C(T ) is almost size independent. Therefore, it is expected that a highesttemperature peak at T ∼ 0.5 shown in Fig. 5 hardly changes even in the thermodynamic limit.
At θ = 0 and θ = 0.1π, we obtain two or three peaks for T < 0.2 in all size. This is consistent with the previous studies for θ = 0 [46][47][48][49]51 . These low-temperature peaks are strongly size-dependent, and thus C(T ) in the thermodynamic limit still unresolved.
At θ = 0.2π, C(T ) exhibits a clear double peak, which has hardly difference between N = 30 and N = 36. Therefore, the existence of this double peak is strongly expected even in the thermodynamic limit at θ = 0.2π.
In addition, the entropy shows a tendency toward a plateau around S(T )/N ∼ 0.3 ∼ ln (2)  However, the origin of the plateau is different, which will be discussed in Sec.V.
At θ = 0.5π (Kitaev limit), S(T ) for N = 30 and N = 36 becomes finite at the lowest temperature (T = 0.0001), being consistent with two-fold (four-fold) degeneracy in the ground state for N = 30 (36). This degeneracy is partially consistent with a previous result using the cluster mean-field method 63 , predicting 2 3L -fold degeneracy in the thermodynamic limit (L → ∞), where L is the linear system size giving the total lattice sites N = 3 × L 2 .
To explore the origin of the multiple-peak structure in S(T ), we calculate S z q (T ) for N = 36 by using the RFTLM, and the results are shown in Fig. 6. When S z q (T ) has the largest intensity at the corner (the edge center) of the extended first Brillouin zone, a √ 3 × √ 3 state (a q = 0 state) appears with short-range order (SRO). At θ = 0 (Heisenberg limit), we obtain a crossover from the paramagnetic state to the √ 3 × √ 3 SRO, and to the q = 0 SRO state, from high to low temperatures. This is the same result obtained by Shimokawa and Kawamura by using the Hams-de Raedt method 49 . At θ > 0 and T = 0.5 where the high-temperature peak of S(T ) appears, we can see that S z q (T ) has zigzag or linear distribution in intensity along the q y direction on q x /π = ±2. This result indicates that the origin of the high-temperature peak is attributed to a crossover from the paramagnetic state to the SRO state with zigzag-or linear-intensity distribution on the SSSF. At 0.1π ≤ θ < 0.5π and T ≤ 0.05, S z q (T ) has the strongest intensity at the edge centers. Therefore, we expect that one of the lower-temperature peaks in C(T ) is a signature of the q = 0, 120 • order. At θ = 0.5π (Kitaev limit), intensity distribution of S z q (T ) has a perfect linear structure. This structure has been obtained in classical spin systems using the Monte Carlo method 63 . Therefore, we can conclude that "order-by-disorder" phenomenon does not occur even in the existence of both the quantum and thermal fluctuations.

C. triangular lattice
The classical ground states in the TL are predicted to be the Z 2 vortex crystal state and the nematic state in 0 ≤ θ ≤ 0.5π [65][66][67][68][69] . We perform finite-temperature calculations for the quantum triangular system. In a recent study, it has been predicted that the specific heat at θ = 0 (Heisenberg limit) has two anomalies at T ∼ 0.2 and T ∼ 0.55 26 . In our calculated C(T ) at θ = 0, a clear peak is obtained at T ∼ 0.2, and a shoulder-like anomaly is obtained at T ∼ 0.6, shown in Fig. 7. A good agreement with the previous work corroborates the validity of our method. In addition, we obtain a gradual change from shoulder-like anomaly to a peak as θ is increased keeping the temperature unchanged. On the other hand, at T ∼ 0.2 and θ = 0.5π for N = 36, the low-temperature peak structure is suppressed. Since this peak exhibits a large-size effect, the specific heat at low temperatures in the thermodynamic limit still remain as a unresolved problem. The entropy of the triangular system is different from that of the kagome system, because there is no plateaulike anomaly in any θ and all N . When θ ≥ 0.375π and N = 36, the ground state has two-fold degeneracy. For this reason, the S(T )/N converges to a value of ln(2)/36 at the lowest temperature T = 0.001 as shown in Fig. 7.
We calculate S z q (T ) of the triangular system for the N = 36 cluster that has a good rotational symmetry  as shown in Fig. 1. Similar to the kagome system, for θ ≥ 0.25π the intensity distribution of S z q (T ) exhibits a zigzag or linear structure along the q y axis on q x /π = ±1 at T = 0.5 where the high-temperature peak in C(T ) appears. For this reason, the high-temperature peak of C(T ) is expected to be a signature of a crossover from the paramagnetic state to a SRO state having the zigzag or linear structure as in the kagome system.
Next we focus on S z q (T ) at T = 0. At θ = 0, S z q (0) has maximum intensity at the corners of the Brillouin zone, which corresponds to the 120 • order as shown in Fig. 9(a). The existence of the 120 • order is consistent with other studies. At θ = 0.5π, S z q (0) has maximum intensity at q = (π, π/ √ 3) and q = (π, −π/ √ 3), meaning the x-stripy order and y-stripy order, respectively, as shown in Fig. 9(b). In the classical system, the ground state has linear-intensity distribution in the SSSF, which is nematic 66 . Therefore, we believe that the "order-bydisorder" phenomenon occurs in the S = 1/2 TL Kitaev model due to the quantum fluctuation. This order has been predicted in the analysis by the linked-cluster expansion and spin-wave theory 86 .

V. DISCUSSION
We compare the results of the kagome and triangular KH model with the honeycomb Kitaev model. In the honeycomb Kitaev model, it has been elucidated that the specific heat has a double-peak structure. In the kagome and triangular KH model, we have found the multiplepeak structures in this work. However, the origin of the double-peak and multiple-peak structures is different. In the honeycomb Kitaev model, the double peak is caused by the itinerant Majorana fermions and Z 2 fluxes freezing at different temperatures 58 . In the kagome system at θ > 0, the high-temperature peak is a consequence of a crossover from the paramagnetic state to a SRO state whose SSSF has a zigzag-or linear-intensity distribution, and one of low-temperature peaks has been expected to be a signature of the q = 0, 120 • order except for θ = 0.5π. In the triangular system, the double-peak structure at θ > 0.25π has the same origin as the kagome system.
The kagome and triangular systems have a significant difference at θ = 0.5π (Kitaev limit). In the triangular system, "order-by-disorder" due to the quantum fluctuations occurs in common with many frustrated quantum spin systems, and the ground state becomes the stripe order. On the other hand, it does not occur in the kagome system.
We have developed new improved FTLMs that are RFTLM and OFTLM. These FTLMs improve the accuracy for all physical quantities at low temperatures compared to the standard FTLM.

VI. SUMMARY
Inspired by the remarkable development of the quantum Kitaev-Heisenberg models in recent years, we investigated the finite-temperature properties of the S = 1/2 KH models on the kagome lattice and triangular lattice by means of improved finite-temperature Lanczos methods. We obtained the multiple peaks in the specific heat in both lattice models. The origin of the hightemperature peak of the specific heat is attributed to a crossover from the paramagnetic state to the SRO state with zigzag or linear structure on the SSSF. We also reveal that at θ = 0.5π (Kitaev limit) in the triangular system, the "order-by-disorder" phenomenon due to the quantum fluctuations occurs, and the ground state exhibits the stripe order. On the other hand, in the kagome system it does not occur even in the presence of both the temperature and quantum fluctuations. We believe this effect is peculiar to the Kitaev model on the kagome lattice.
We have succeeded in improving the finite-temperature Lanczos method. For larger systems, we can expect further improvements, especially speed-up calculations, using a technique for decomposing full Hilbert space with several symmetries such as the case of SPINPACK 87 . Next target for finite-temperature calculations will be lattices with 48 sites, which remains as a future work.

ACKNOWLEDGMENTS
This work was supported by MEXT, Japan, as a social and scientific priority issue (creation of new functional devices and high-performance materials to support nextgeneration industries) to be tackled by using a post-K computer. The numerical calculation was carried out at the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo and at the Yukawa Institute Computer Facility, Kyoto University.