# Ground-state and thermodynamic properties of an Kitaev model

\abstWe study ground-state and thermodynamic properties of an Kitaev model. We first clarify the existence of global parity symmetry in addition to the local symmetry on each plaquette, which enables us to perform the large scale calculations up to 24 sites. It is found that the ground state should be singlet and its energy is estimated as , where is the Kitaev exchange coupling. We find that a lowest excited state belongs to the same subspace as the ground state and the gap monotonically decreases with increasing system size, which suggests that the ground state of the Kitaev model is gapless. By means of the thermal pure quantum states, we clarify finite temperature properties characteristic of the Kitaev models with .

Frustrated quantum spin systems have been the subjects of considerable interest. One of the intriguing systems is the quantum spin Kitaev model on the honeycomb lattice [1], where strong frustration emerges due to direction-dependent Ising interactions. It is known that this model is exactly solvable and its ground state is a quantum spin liquid state, where spin degrees of freedom are decoupled by itinerant Majorana fermions and fluxes [1, 2, 3, 4]. Two energy scales for distinct degrees of freedom yield interesting finite temperature properties such as double peak structure in the specific heat, spin dynamics, and thermal Hall coefficients at low temperatures [5, 6, 7, 8, 9, 10]. Since magnetic properties of certain Mott insulators with strong spin-orbit coupling should be described by the Kitaev model and its extentions [11], a lot of theoretical studies for these models [12, 13, 14, 15, 16, 17] and experimental studies on candidates materials such as IrO(=Na, Li) [18, 19, 20, 21], -RuCl [22, 23, 24, 25] and [26] have been done so far.

The key of these features characteristic of the Kitaev model is the existence of the local conserved quantity defined at each plaquette in the honeycomb lattice, which is responsible for the absence of the long-range spin correlations. From a fundamental viewpoint of quantum spin systems, a question arises whether or not such interesting properties appear in generalized Kitaev models with arbitrary spins. It is known that in the general spin Kitaev model [27], there exist local conserved quantities and the ground state is nonmagnetic. In spite of the presence of the local invariants, details of its excitations and finite temperature properties remain unclear. Very recently thermodynamic properties in the Kitaev model with spin has been examined for an 8-site cluster [17]. It has been found that, in the specific heat, double peak structure appears only for the and models, while a single peak appears in the model. However, the system size may not be large enough to discuss thermodynamic properties for the Kitaev models on the honeycomb lattice. In addition, the role of the local conserved quantities was still unclear. Therefore, it is necessary to clarify how the degrees of freedom inherent in the Kitaev model affects ground state and finite temperature properties.

In this Letter, we mainly focus on the ground state and finite temperature properties of the Kitaev model. First, we show the existence of the global parity symmetry in addition to the symmetry in each plaquette in the system [27]. Using these conserved quantities, we examine ground state and finite-temperature properties in the large clusters up to 24 sites. We clarify that ground state has no degeneracy and belongs to the subspace with zero-fluxes, which is consistent with the semiclassical results [27]. We find that a lowest excited state belongs to the same subspace as the ground state and the gap monotonically decreases with increasing system size, which suggests that the ground state of the Kitaev model is gapless. Furthermore, using the thermal pure quantum (TPQ) state methods, we demonstrate that the multiple peak structure appears in the specific heat, similar to the Kitaev model [5]. Thermodynamic properties in the generic spin- Kitaev models with are also addressed.

We consider the Kitaev model given by

(1) |

where is an component of an spin operator at the th site in the honeycomb lattice. The ferromagnetic interaction is defined on three different types of the nearest neighbor bonds, (red), (blue), and bonds (green), respectively [see Fig. 1(a)].

In this model, the basis sets are convenient [28], which is explicitly given as

(2) | ||||

(3) | ||||

(4) |

where is an eigenstate of the -component of the spin operator with an eigenvalue . In this representation, we obtain the following relations as

(5) |

Now, we consider the local Hamiltonian for the bond between and sites, . When a certain state includes an state in th and/or th sites, due to eq. (5). Then, nonzero components in the Hamiltonian are explicitly given as

(6) | ||||

(7) | ||||

(8) | ||||

(9) |

We here note that the numbers of and states are not conserved, but their parity is conserved. When one considers all Kitaev interactions in the system, there exist global parity symmetries for the number of states. In other word, the Kitaev Hamiltonian eq. (1) commutes the parity operators for states as

(10) |

where is the number operator for states, and , and commute with each other. In addition, there exists a local conserved quantity defined at each plaquette composed of the sites labeled as [see Fig. 1(b)] [27]:

(11) |

This operator commutes with not only on another plaquette and the Hamiltonian eq. (1), but also the parity operators . Then, the Hilbert space of the Hamiltonian can be classified into each subspace specified by , where and are eigenvalues of the operators and . This enables us to perform the numerical calculations for large clusters in the Lanczos and TPQ state methods. For example, the Hilbert space in the cluster with is divided into 2,048 subspaces, in which the maximum rank of the block-diagonalized Hamiltonian is 140,279,025. Applying the Lanczos and TPQ state methods to the systems with periodic boundary conditions, we discuss ground state and finite temperature properties.

In our study, we deal with the Kitaev model on the clusters with , and 24, which are shown in Fig. 1(a). Calculating the lowest energy for each cluster, we find that the ground state has no degeneracy and always belongs to the subspace , corresponding to the zero-flux sector. This is consistent with the results obtained from the semiclassical approach [27]. The ground state energy for each cluster is shown in Fig. 2(a).

In the thermodynamic limit, the ground state energy appears to be given by although the systematic finite size scaling is necessary to be deduced quantitatively.

We also discuss the elementary excitations in the Kitaev model, examining all subspaces for the clusters. Figure 2(b) shows lowest excitation energies in the subspaces and , where is the subspace with eigenvalues , and means the complement space of . The excited states strongly depend on the size and/or shape of clusters, as shown in Fig. 2(b). Therefore, we believe that the largest cluster with isotropic geometry should capture the nature of the excitations in the thermodynamic limit. In the system, the lowest excited state belongs to the subspace , which is the same subspace for the ground state. Furthermore, we find that the excitation gap in the subspace monotonically decreases with increasing the system size and is much smaller than the exchange constant . These suggest that in the thermodynamic limit, the Kitaev model has gapless excitations in the same subspace . This is similar to the case with the Kitaev model, where the gapless dispersion for Majorana fermions is obtained in the subspace corresponding to the zero-flux sector [1].

Next, we examine thermodynamic properties in the Kitaev model. We calculate thermodynamic quantities for a small cluster with , obtaining all eigenenergies by means of the exact diagonalization method. Furthermore, we make use of the TPQ states for larger clusters. As for two kinds of clusters with , we evaluate thermodynamic quantities and expectation values by means of hundred independent TPQ states. By contrast, the numerical cost is expensive for the site cluster, and only two TPQ states are treated. The obtained results are shown in Fig. 3.

We clearly find in Fig. 3(a), the multiple peak structure in the specific heat , which is consistent with the previous work [17]. Now, we focus on the peak structure around . Since the specific heat around is almost independent of the cluster size, this peak structure appears in the thermodynamic limit. On the other hand, low temperature structure is sensitive to the cluster size, which originates from the fact that low-energy excitations depend on the size and/or shape of clusters [see Fig. 2(b)]. Nevertheless, the peak appears around except for the smallest cluster with , which suggests that the low temperature peak indeed appears in the thermodynamic limit. In addition to these features, the shoulder or peak structure appears in lower energy region , as shown in Fig. 3(a). Since it strongly depends on the clusters, further detailed analysis with larger clusters is necessary to clarify whether or not such low-temperature behavior appears in the thermodynamic limit.

We also calculate the entropy , and the expectation values of the local conserved quantity , and spin-spin correlations , which are defined as

(12) | ||||

(13) |

where is the number of plaquettes, and is the number of bonds. Note that is proportional to the energy as . The results are shown in Figs. 3(b) and (c). It is found that, decreasing temperatures, nearest-neighbor spin-spin correlations develop around and the plateau behavior appears in the curve of the entropy at around . On the other hand, the expectation value changes rapidly at lower temperature . Therefore, the peaks in the specific heat at and should be related to spin-spin correlations and local conserved quantity.

From the above discussions, we have confirmed that, in the Kitaev model, at least, double peak structure appears in the specific heat. This originates from two distinct temperature dependences for nearest-neighbor spin-spin correlations and the local conserved quantity, which are essentially the same as those in the Kitaev model [5]. Then, similar finite-temperature behavior is naively expected even in larger spin models since they have the local conserved quantity [27]. To clarify this, we also examine finite-temperature properties for the Kitaev models with and , using the TPQ states for the clusters with and , respectively. First, we have compared the results for the and clusters in the spin system, and have confirmed that the higher temperature peak is not changed in the specific heat. Therefore, we believe that the site calculations can capture higher temperature peak in the thermodynamic limit. Meanwhile, the characteristic temperature for low temperature peak in the specific heat depends on the system size, which is similar to the case (see Fig. 3).

Figure 4 shows the specific heat and entropy in the Kitaev models with . Namely, the results are obtained by the Monte Carlo method for sites [5, 8]. The results for and are obtained by the TPQ state methods for , and clusters, respectively. We find multiple peak structures in the specific heat for the spin- Kitaev model in Fig. 4(a). This is contrast to the previous work [17], where a single peak is observed in the specific heat for the cases in an 8-site cluster. This should originate from the fact that the system size treated in previous works is too small to discuss finite temperature properties in the thermodynamic limit. Now, we focus on the higher-temperature peak in the specific heat, which can be described quantitatively in our methods. It is found that the characteristic temperature and the height of the peak monotonically increase with increasing the amplitude of spin , as shown in Fig. 4(a). To clarify the spin dependence, we also plot its characteristic temperature and the value of the specific heat at in the inset of Fig. 4(a). We numerically find that these quantities are almost scaled as and . These are in contrast to the spin dependence of the characteristic temperature derived from the high-temperature expansion: the specific heat is expanded as with . The difference in the characteristic temperatures implies that the fractionalization of spin occurs in the spin- Kitaev models as well as the Kitaev model [1]. To confirm this, we calculate the temperature dependence of the entropy, as shown in Fig. 4(b). We find that the plateau appears in the curve of the entropy at , suggesting that spin with the amplitude is fractionalized into two halves with distinct energy scales. This is corroborated by the fact that the temperature of normalized entropy is well scaled by , as shown in the inset of Fig. 4(b). We naively expect that the decrease of temperatures induces another peak in the specific heat, where the residual entropy is released. Unfortunately, the low-temperature peak could not be discussed quantitatively since the clusters treated in the present study are not enough large to describe low energy excitations correctly. Therefore, the systematic analysis with larger clusters are needed to clarify whether peak and/or shoulder structures appears or not in the lower temperature region, which is left for a future work.

In summary, we have studied an Kitaev model to discuss ground-state and finite temperature properties. We have demonstrated the existence of parity symmetry as well as the local conserved quantity, which enables us to perform the calculations in large clusters up to twenty-four sites. We have found that the ground state is singlet and its energy is estimated as . In addition, a lowest excited state belongs to the same subspace as the ground state and the gap monotonically decreases with increasing system size, which suggests that the ground state of the Kitaev model is gapless. By means of the TPQ states, we have examined finite temperature properties and clarified that the multiple peak structure in the specific heat appears in the Kitaev models with spin . This should be a common feature of the general Kitaev model, which might originate from the presence of a local conserved quantity on each plaquette in the lattice.

This work is supported by Grant-in-Aid for Scientific Research from JSPS, KAKENHI Grant Number JP17K05536, JP16H01066 (A.K.) and JP16K17747, JP16H02206, JP16H00987 (J.N.). Parts of the numerical calculations were performed in the supercomputing systems in ISSP, the University of Tokyo.

## References

- [1] A. Kitaev: Ann. Phys. (N. Y.) 321 (2006) 2.
- [2] X.-Y. Feng, G.-M. Zhang, and T. Xiang: Phys. Rev. Lett. 98 (2007) 087204.
- [3] H.-D. Chen and J. Hu: Phys. Rev. B 76 (2007) 193101.
- [4] H.-D. Chen and Z. Nussinov: J. Phys. A 41 (2008) 075001.
- [5] J. Nasu, M. Udagawa, and Y. Motome: Phys. Rev. B 92 (2015) 115122.
- [6] Y. Yamaji, T. Suzuki, T. Yamada, S.-i. Suga, N. Kawashima, and M. Imada: Phys. Rev. B 93 (2016) 174425.
- [7] J. Yoshitake, J. Nasu, and Y. Motome: Phys. Rev. Lett. 117 (2016) 157203.
- [8] J. Nasu, J. Yoshitake, and Y. Motome: Phys. Rev. Lett. 119 (2017) 127204.
- [9] T. Suzuki and S.-i. Suga: (unpublished) arXiv:1802.00545.
- [10] Y. Yamaji, T. Suzuki, and M. Kawamura: (unpublished) arXiv:1802.02854.
- [11] G. Jackeli and G. Khaliullin: Phys. Rev. Lett. 102 (2009) 017205.
- [12] J. c. v. Chaloupka, G. Jackeli, and G. Khaliullin: Phys. Rev. Lett. 110 (2013) 097204.
- [13] Y. Yamaji, Y. Nomura, M. Kurita, R. Arita, and M. Imada: Phys. Rev. Lett. 113 (2014) 107201.
- [14] T. Suzuki, T. Yamada, Y. Yamaji, and S.-i. Suga: Phys. Rev. B 92 (2015) 184411.
- [15] R. R. P. Singh and J. Oitmaa: Phys. Rev. B 96 (2017) 144414.
- [16] A. Koga, S. Nakauchi, and J. Nasu: preprint (unpublished) arXiv:1705.09659.
- [17] T. Suzuki and Y. Yamaji: Physica B (in press.).
- [18] Y. Singh and P. Gegenwart: Phys. Rev. B 82 (2010) 064412.
- [19] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart: Phys. Rev. Lett. 108 (2012) 127203.
- [20] R. Comin, G. Levy, B. Ludbrook, Z.-H. Zhu, C. N. Veenstra, J. A. Rosen, Y. Singh, P. Gegenwart, D. Stricker, J. N. Hancock, D. van der Marel, I. S. Elfimov, and A. Damascelli: Phys. Rev. Lett. 109 (2012) 266406.
- [21] S. K. Choi, R. Coldea, A. N. Kolmogorov, T. Lancaster, I. I. Mazin, S. J. Blundell, P. G. Radaelli, Y. Singh, P. Gegenwart, K. R. Choi, S.-W. Cheong, P. J. Baker, C. Stock, and J. Taylor: Phys. Rev. Lett. 108 (2012) 127204.
- [22] K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim: Phys. Rev. B 90 (2014) 041112.
- [23] Y. Kubota, H. Tanaka, T. Ono, Y. Narumi, and K. Kindo: Phys. Rev. B 91 (2015) 094422.
- [24] J. A. Sears, M. Songvilay, K. W. Plumb, J. P. Clancy, Y. Qiu, Y. Zhao, D. Parshall, and Y.-J. Kim: Phys. Rev. B 91 (2015) 144420.
- [25] M. Majumder, M. Schmidt, H. Rosner, A. A. Tsirlin, H. Yasuoka, and M. Baenitz: Phys. Rev. B 91 (2015) 180401.
- [26] S. Bette, T. Takayama, K. Kitagawa, R. Takano, H. Takagi, and R. E. Dinnebier: Dalton Trans. 46 (2017) 15216.
- [27] G. Baskaran, D. Sen, and R. Shankar: Phys. Rev. B 78 (2008) 115116.
- [28] H. Tomishige, J. Nasu, and A. Koga: Phys. Rev. B 97 (2018) 094403.