Finite-temperature phase transition to a Kitaev spin liquid phase on a hyperoctagon lattice: A large-scale quantum Monte Carlo study
Abstract
The quantum spin liquid is an enigmatic quantum state in insulating magnets, in which conventional long-range order is suppressed by strong quantum fluctuations. Recently, an unconventional phase transition was reported between the low-temperature quantum spin liquid and the high-temperature paramagnet in the Kitaev model on a three-dimensional hyperhoneycomb lattice. Here, we show that a similar “liquid-gas” transition takes place in another three-dimensional lattice, the hyperoctagon lattice. We investigate the critical phenomena by adopting the Green-function based Monte Carlo technique with the kernel polynomial method, which enables systematic analysis of up to sites. The critical temperature is lower than that in the hyperhoneycomb case, reflecting the smaller flux gap. We also discuss the transition on the basis of an effective model in the anisotropic limit.
pacs:
Valid PACS appear hereI Introduction
The quantum spin liquid (QSL) is an exotic magnetic state of matter in solids, in which long-range ordering is suppressed by strong quantum fluctuations even in the ground state Anderson (1973); Fazekas and Anderson (1974); Balents (2010). As a consequence of intensive study over several decades, substantial progress has been recently made in experimental search for QSLs Kanoda and Kato (2011); Gingras and McClarty (2014); Okamoto et al. (2007); Coldea et al. (2003); Hiroi et al. (2001); Shores et al. (2005); Okamoto et al. (2009); Han et al. (2014). In such experiments, however, QSLs are usually identified by the absence of any phase transition to a long-range ordered state. In other words, it is supposed that the low-temperature () QSL is adiabatically connected to the high- paramagnetic state. This common belief relies on the magnetic analog of liquid helium which remains as a liquid down to the lowest due to strong quantum fluctuations Anderson (1973). Nonetheless, in conventional fluids, gas and liquid can be distinguished by a phase transition. Hence, a similar phase transition can also be expected in the case of magnetic states of matter, namely, between the “spin gas” (paramagnet) and “spin liquid” (QSL). However, the thermodynamics of QSLs remains elusive, despite the crucial importance for understanding of existing and forthcoming experiments.
Recently, an exotic finite- phase transition between the low- QSL and the high- paramagnet was reported for an extension of the Kitaev model defined on a three-dimensional (3D) hyperhoneycomb lattice as well as its anisotropic limit Nasu et al. (2014a, b). The Kitaev model is an exactly-soluble quantum spin model, whose ground state provides a canonical example of QSLs Kitaev (2006), and has attracted attention as it might describe the magnetism in some spin-orbit entangled Mott insulatorsJackeli and Khaliullin (2009). The finding of the exotic phase transition was brought by a newly-developed quantum Monte Carlo (QMC) method based on a Majorana fermion representation, which is free from the negative sign problem. This phase transition is not accompanied by any symmetry breaking, but can be explained by a proliferation of loops for excited fluxes Kimchi et al. (2014); Nasu et al. (2014b, a). The emergent loop degree of freedom is specific to 3D; the finite- phase transition is absent and turns into a crossover in two dimensions Nasu et al. (2015). The surprising result urges reconsideration of both experimental and theoretical quests for QSLs.
As demonstrated in the hyperhoneycomb case Mandal and Surendran (2009), the Kitaev model can be extended to any tri-coordinate lattices Yang et al. (2007); O’Brien et al. (2016); Hermanns and Trebst (2014); Yao and Kivelson (2007). Although all the extensions retain the QSL nature Baskaran et al. (2007), the exact ground state (the spatial configuration of fluxes) is obtained for some limited cases because Lieb’s theorem is not applicable to generic tri-coordinate lattices Lieb (1994). Moreover, finite- properties have not been studied in most cases. As there are a variety of 3D tri-coordinate lattices O’Brien et al. (2016); Wells (1977), it is interesting to investigate such 3D models for clarifying the universal aspects of the phase transition between the QSL and the paramagnet and for exploring further exotic phase transitions.
In this paper, we investigate finite- properties in the Kitaev model on another 3D lattice, the hyperoctagon lattice Hermanns and Trebst (2014) (see Fig. 1). For this model, the exact solution is not available, as Lieb’s theorem cannot be applied. Improving the QMC method by employing the Green function technique with the kernel polynomial method (KPM), we compute much larger systems compared to the previous studies. From the systematic analysis of clusters up to sites, we find that the hyperoctagon model exhibits another example of finite- phase transitions between the QSL and the paramagnet. We find that the critical temperature estimated from the careful analysis of the finite-size effects is lower than that for the hyperhoneycomb case, reflecting the smaller flux gap O’Brien et al. (2016); Kimchi et al. (2014). In addition, we show that the anisotropic limit of the hyperoctagon model becomes equivalent to that of the hyperhoneycomb one, which supports the common origin of the phase transition on two lattices.
The paper is structured as follows. In Sec. II, we introduce the Kitaev model on a hyperoctagon lattice, and briefly describe the numerical method used in the previous studies to analyze the finite- properties of the Kitaev model. Then we introduce the improved QMC method by using the Green-function based KPM (GFb-KPM) and present the benchmark results. In Sec. III, we present numerical results. We identify the finite- phase transition and estimate the critical temperature from careful analyses of finite-size effects for three different physical quantities. In Sec. IV, we discuss the origin of the phase transition on the hyperoctagon lattice, with a consideration of the anisotropic limit of the model in comparison with that for the hyperhoneycomb case. We also discuss a correlation between the critical temperature and the flux gap, with a brief comment on the ground state. Finally, Sec. V is devoted to the summary.
Ii Model and Method
We study an extension of the Kitaev model Kitaev (2006) to a 3D hyperoctagon lattice shown in Fig. 1 Hermanns and Trebst (2014). The Hamiltonian is given by
(1) |
where , , and are the Pauli matrices describing a spin- state at site . The sum is taken over the nearest-neighbor (NN) sites on the three different types of bonds, which are shown by red (), green (), and blue () in Fig. 1; is the exchange constant for each bond. Similar to the original Kitaev model on a honeycomb lattice Kitaev (2006), the hyperoctagon model has a conserved quantity called the fluxes for each ten-site plaquette (see Fig. 1), , where the product is taken for all the bonds comprising the plaquette. While the ground state of this model is not exactly obtained, it was deduced to be a flux-free state, where all are Hermanns and Trebst (2014), similar to the honeycomb Kitaev (2006) and hyperhoneycomb cases Mandal and Surendran (2009).
We investigate thermodynamic properties of the model in Eq. (1) by the QMC simulation based on a Majorana fermion representation Nasu et al. (2014a). In this technique, the lattice is regarded as an assembly of one-dimensional chains composed of two types of bonds and the Jordan-Wigner transformation is applied to each chain Chen and Hu (2007); Feng et al. (2007); Chen and Nussinov (2008). By taking such chains along the direction composed of the red and green bonds as in Fig. 1, the Kitaev Hamiltonian in Eq. (1) is rewritten by two types of Majorana fermions and as
(2) |
where is a local conserved quantity taking ( is the index for the bond); it satisfies , where the product is taken for all the bonds included in the ten-site loop. The sum of is taken for the NN sites on the bonds with () if the chain begins with () (see Fig. 1). From the form of Eq. (2), the system is regarded as free Majorana fermions coupled with the variables . Hence, one can compute thermodynamic properties by Monte Carlo (MC) sampling over the configurations with the statistical weight , where is the free energy of the Majorana fermions for a configuration , and is the inverse temperature (we set the Boltzmann constant ). As is positive definite, the QMC simulation is free from the negative sign problem Nasu et al. (2014a); Nasu and Motome (2015); Nasu et al. (2015, 2016).
In the previous studies Nasu et al. (2014a); Nasu and Motome (2015); Nasu et al. (2015, 2016), was calculated by the exact diagonalization (ED), which leads to the computational cost of for one MC sweep where is the number of spins. This has limited the reachable system size to less than in a reasonable computational time. In order to reduce the computational cost and enable systematic analysis of up to larger , in the present study, we apply the GFb-KPM to estimate Weiße (2009). In this method, a change of the Majorana fermion density of states in a MC update is directly computed from a small number of Green functions obtained by the KPM Weiße et al. (2006), which reduces the computational cost to .
Let us describe the method briefly. Suppose if the one-body Majorana Hamiltonian for a given configuration , , is modified to by a local flip of , then the energy spectrum is given by the solution of , where is the Green function satisfying ( is the unit matrix). By extending to a complex function by , the difference in the free energy by the local flip is given by
(3) |
can be calculated in a compact manner by using a few components of Green’s functions Weiße et al. (2006). The Green functions are obtained by using the KPM with Chebyshev polynomials. For instance, the onsite component is obtained as
(4) |
where is the th diagonal element of the th Chebyshev moment of rescaled by a factor of to fit the eigenvalues in the range of , and is the cutoff of the expansion Weiße (2009). is calculated efficiently by sparse matrix-vector multiplications Weiße et al. (2006).
In the following calculations, we focus on the isotropic case and set the energy unit as . We consider the clusters from ( unit cells) to sites ( unit cells) with the open boundary condition in the direction and the periodic ones in the remaining two directions. We typically perform ten independent MC runs with MC steps for measurements after to steps for thermalization. In addition to MC sweeps by a single flip of , we adopt the replica exchange MC technique Hukushima and Nemoto (1996).
Figure 2 shows the benchmark of the GFb-KPM method. We compare the QMC data of the specific heat per site, , computed
by the ED method and the GFb-KPM method for the -site cluster. For the GFb-KPM calculations, we present the data with a
different cutoff of the Chebyshev expansion, , from to . Figure 2(a) shows dependence of . We find
that the GFb-KPM data well converge to the ED ones, even for the smallest , except for the low- region around the low- peak
of . Figure 2(b) shows the convergence of the GFb-KPM data in terms of in the low- region: the symbols indicate
the GFb-KPM data for each , while the horizontal lines the ED ones (the dashed lines indicate the statistical errorbars). We find that
the GFb-KPM results quickly converge to the ED ones while increasing . On the basis of this benchmark as well as the fact that the
GFb-KPM converges faster for larger system sizes Weiße et al. (2006), we take for the -site cluster and
for larger clusters in the following GFb-KPM calculations
Iii Results
Figure 3(a) shows the QMC results of the specific heat per site, . Similar to the hyperhoneycomb case Nasu et al. (2014a), exhibits two peaks. The higher- peak is closely related with the development of NN spin correlations shown in Fig. 3(c), which corresponds to the lowering of the kinetic energy of the itinerant Majorana fermions Nasu et al. (2015, 2014a). On the other hand, the lower- peak in is associated with the coherent alignment of the variables Nasu et al. (2015, 2014a); see the thermal average plotted in Fig. 3(c) ( is the number of ten-site plaquettes in the system). Figure 3(b) shows the results of the entropy per site, . At each peak of , a half of entropy is released; the higher- release comes from itinerant Majorana fermions and the lower one comes from variables . Similar behavior was reported for the Kitaev models on the honeycomb and hyperhoneycomb lattices Nasu et al. (2015, 2014a). This is the common feature in the Kitaev QSLs called thermal fractionalization of quantum spins Nasu et al. (2015).
In Fig. 3(a), while the high- peak of is almost system-size independent, the low- peak exhibits significant system-size dependence: the peak becomes sharper and shifts to lower while increasing the system size. The enlarged view of the low- part is shown in Fig. 4(a). The system-size dependence is indicative of a phase transition, as in the hyperhoneycomb case Nasu et al. (2014a). We note that the height of the peaks does not show systematic behavior, which is presumably due to the different shapes of clusters. We find similar critical behavior in thermal fluctuations of , defined by
(5) |
As shown in Fig. 4(b), shows a similar peak to . We note that this quantity gives a measure of the specific heat in the case of the anisotropic limit of the Kitaev model where the Hamiltonian is described by only; namely, gives a measure of energy fluctuations related to . Thus, the common critical behavior between and indicates that the phase transition is driven by the variables .
In order to characterize the phase transition, following the previous study for the hyperhoneycomb case Nasu et al. (2014a), we also
measure the Wilson loop on a plane perpendicular to one of the periodic boundaries
From the critical behaviors of , , and , we estimate the critical temperature . The insets in Fig. 4 present the system-size extrapolation of the peak of and and the inflection point of dependence of . The extrapolated values to well coincide with each other: by , by , and by (the numbers in the parentheses represent the errors in the last digits). This good agreement indicates the validity of our numerical simulation and the system-size analysis.
Iv Discussion
From the comparison to the previous study for the hyperhoneycomb model Nasu et al. (2014a), we conclude that the phase transition on
the hyperoctagon lattice is another realization of the “liquid-gas” phase transition between the low- QSL and the high- paramagnet.
Similar to the hyperhoneycomb case, we do not observe any indication of discontinuity in the transition
The common phase transition between the hyperoctagon and hyperhoneycomb cases is also inferred by considering the anisotropic limit of the models, and . For both cases, the low-energy effective model is obtained by the similar perturbation theory to that used for the derivation of the toric code from the Kitaev model in Ref. Kitaev, 2006. In the perturbation theory, the original Hilbert space is projected onto the low-energy subspace where each dimer made of two spins connected by the strong bond, say and , takes only two states out of four; these two states are represented by a pseudospin . When one defines the pseudospins at the centers of bonds, the lattice is reduced to the diamond lattice for both hyperoctagon and hyperhoneycomb cases. For the hyperhoneycomb case, the low-energy effective model was obtained in Ref. Mandal and Surendran, 2009 by considering the eighth-order perturbation. The effective model is described by the projected fluxes that are defined for each six-site plaquette in the diamond lattice, and has the form of , where is either or depending on the number of bond included in the plaquette . It was numerically demonstrated that the effective model exhibits a finite- phase transition between the high- paramagnet and the low- QSL Nasu et al. (2014b). Performing the similar procedure, we find that the low-energy effective model for the hyperoctagon case has the same form, except for the sign of the coupling constant . As the opposite sign does not cause any difference in the thermodynamics, the effective model shows the same finite- phase transition as in the hyperhoneycomb case. Although this correspondence is only in the anisotropic limit, it strongly suggests a common mechanism of the finite- phase transition between the two models even in the isotropic case.
On the other hand, we note that our estimate of for the hyperoctagon model is lower than that for the hyperhoneycomb case, as shown in Table 1 Nasu et al. (2014a). As the phase transition is caused by the loop proliferation, is closely related with the loop tension Kimchi et al. (2014). The loop tension is determined by the flux gap , which is defined by the minimum excitation energy for flipping from the ground state. The flux gap is estimated as for the hyperoctagon lattice and for the hyperhoneycomb lattice O’Brien et al. (2016). The values of and summarized in Table 1 indicate a good correlation between and .
Finally, let us briefly comment on the ground state. As shown in Fig. 3(c), our QMC data for converge to at low . This suggests that the system is likely to be in the flux-free state with all when , as deduced by variational arguments in the previous study Hermanns and Trebst (2014).
Lattice | Critical temperature | Flux gap |
---|---|---|
Hyperoctagon | 0.00405(10) | 0.030(3) |
Hyperhoneycomb | 0.00519(9) | 0.043(3) |
V Summary
In summary, we have investigated the “liquid-gas” phase transition from the low- QSL to the high- paramagnet on the extension of the Kitaev model to the 3D hyperoctagon lattice. Using an improved QMC method by the Green function technique with the KPM, we successfully dealt with much larger-size systems compared to the previous studies. From the systematic analysis of clusters up to sites, we found a finite- phase transition similar to the previous hyperhoneycomb case. The result was also supported by considering the anisotropic limit of the model on two lattices. We also showed that the critical temperature on the hyperoctagon lattice is lower than in the hyperhoneycomb case, reflecting the smaller flux gap in the former case. It will be interesting to extend the present study to other 3D Kitaev models for an exploration of further exotic phase transitions.
Acknowledgements.
The authors thank J. Nasu, R. Ozawa, and J. Yoshitake for fruitful discussions. This research was supported by JSPS Grants-in-Aid for Scientific Research Grants No. JP15K13533, No. JP16H02206, and No. 26800199. Parts of the numerical calculations were performed in the supercomputing systems in ISSP, the University of Tokyo.Footnotes
- preprint: APS/123-QED
- We compute the two smallest clusters with and by QMC with ED.
- We define on a half of the largest loop on the plane which is halved in the periodic boundary condition.
- We carefully examine the histogram of the internal energy, but do not find any clear signature of bifurcation.
References
- P. W. Anderson, “Resonating valence bonds: A new kind of insulator?,” Mater. Res. Bull. 8, 153 (1973).
- P. Fazekas and P. W. Anderson, “On the ground state properties of the anisotropic triangular antiferromagnet,” Phil. Mag. 30, 423 (1974).
- L. Balents, “Spin liquids in frustrated magnets,” Nature (London) 464, 199 (2010).
- K. Kanoda and R. Kato, “Mott Physics in Organic Conductors with Triangular Lattices,” Annu. Rev. Condens. Matter Phys. 2, 167 (2011).
- M. J. P. Gingras and P. A. McClarty, “Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets,” Rep. Prog. Phys 77, 056501 (2014).
- Y. Okamoto, M. Nohara, H. Aruga-Katori, and H. Takagi, “Spin-Liquid State in the Hyperkagome Antiferromagnet NaIrO,” Phys. Rev. Lett. 99, 137207 (2007).
- R. Coldea, D. A. Tennant, and Z. Tylczynski, “Extended scattering continua characteristic of spin fractionalization in the two-dimensional frustrated quantum magnet CsCuCl observed by neutron scattering,” Phys. Rev. B 68, 134424 (2003).
- Z. Hiroi, M. Hanawa, N. Kobayashi, M. Nohara, H. Takagi, Y. Kato, and M. Takigawa, “Spin- Kagome-Like Lattice in Volborthite CuVO(OH)HO,” J. Phys. Soc. Jpn. 70, 3377 (2001).
- M. P. Shores, E. A. Nytko, B. M. Bartlett, and D. G. Nocera, “A Structurally Perfect Kagome Antiferromagnet,” J. Am. Chem. Soc. 127, 13462 (2005).
- Y. Okamoto, H. Yoshida, and Z. Hiroi, “Vesignieite BaCuVO(OH) as a Candidate Spin- Kagome Antiferromagnet,” J. Phys. Soc. Jpn. 78, 033701 (2009).
- T.-H. Han, J. Singleton, and J. A. Schlueter, “Barlowite: A Spin- Antiferromagnet with a Geometrically Perfect Kagome Motif,” Phys. Rev. Lett. 113, 227203 (2014).
- J. Nasu, M. Udagawa, and Y. Motome, “Vaporization of Kitaev Spin Liquids,” Phys. Rev. Lett. 113, 197205 (2014a).
- J. Nasu, T. Kaji, K. Matsuura, M. Udagawa, and Y. Motome, “Finite-temperature phase transition to a quantum spin liquid in a three-dimensional Kitaev model on a hyperhoneycomb lattice,” Phys. Rev. B 89, 115125 (2014b).
- A. Kitaev, “Anyons in an exactly solved model and beyond,” Ann. Phys. (N.Y.) 321, 2 (2006).
- G. Jackeli and G. Khaliullin, “Mott Insulators in the Strong Spin-Orbit Coupling Limit: From Heisenberg to a Quantum Compass and Kitaev Models,” Phys. Rev. Lett. 102, 017205 (2009).
- I. Kimchi, J. G. Analytis, and A. Vishwanath, “Three-dimensional quantum spin liquids in models of harmonic-honeycomb iridates and phase diagram in an infinite- approximation,” Phys. Rev. B 90, 205126 (2014).
- J. Nasu, M. Udagawa, and Y. Motome, “Thermal fractionalization of quantum spins in a Kitaev model: Temperature-linear specific heat and coherent transport of Majorana fermions,” Phys. Rev. B 92, 115122 (2015).
- S. Mandal and N. Surendran, “Exactly solvable Kitaev model in three dimensions,” Phys. Rev. B 79, 024426 (2009).
- S. Yang, D. L. Zhou, and C. P. Sun, “Mosaic spin models with topological order,” Phys. Rev. B 76, 180404(R) (2007).
- K. O’Brien, M. Hermanns, and S. Trebst, “Classification of gapless spin liquids in three-dimensional Kitaev models,” Phys. Rev. B 93, 085101 (2016).
- M. Hermanns and S. Trebst, “Quantum spin liquid with a Majorana Fermi surface on the three-dimensional hyperoctagon lattice,” Phys. Rev. B 89, 235102 (2014).
- H. Yao and S. A. Kivelson, “Exact Chiral Spin Liquid with Non-Abelian Anyons,” Phys. Rev. Lett. 99, 247203 (2007).
- G. Baskaran, S. Mandal, and R. Shankar, “Exact Results for Spin Dynamics and Fractionalization in the Kitaev Model,” Phys. Rev. Lett. 98, 247201 (2007).
- E. H. Lieb, “Flux Phase of the Half-Filled Band,” Phys. Rev. Lett. 73, 2158 (1994).
- A. F. Wells, Three Dimensional Nets and Polyhedra, 1st ed. (John Wiley and Sons Inc., New York, 1977).
- H.-D. Chen and J. Hu, “Exact mapping between classical and topological orders in two-dimensional spin systems,” Phys. Rev. B 76, 193101 (2007).
- X.-Y. Feng, G.-M. Zhang, and T. Xiang, “Topological Characterization of Quantum Phase Transitions in a Spin- Model,” Phys. Rev. Lett. 98, 087204 (2007).
- H.-D. Chen and Z. Nussinov, “Exact results of the Kitaev model on a hexagonal lattice: spin states, string and brane correlators, and anyonic excitations,” J. Phys. A 41, 075001 (2008).
- J. Nasu and Y. Motome, “Thermodynamics of Chiral Spin Liquids with Abelian and Non-Abelian Anyons,” Phys. Rev. Lett. 115, 087203 (2015).
- J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, and R. Moessner, “Fermionic response from fractionalization in an insulating two-dimensional magnet,” Nature Physics 12, 912 (2016).
- A. Weiße, “Green-Function-Based Monte Carlo Method for Classical Fields Coupled to Fermions,” Phys. Rev. Lett. 102, 150604 (2009).
- A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, “The kernel polynomial method,” Rev. Mod. Phys. 78, 275 (2006).
- K. Hukushima and K. Nemoto, “Exchange Monte Carlo Method and Application to Spin Glass Simulations,” J. Phys. Soc. Jpn. 65, 1604 (1996).
- We compute the two smallest clusters with and by QMC with ED.
- We define on a half of the largest loop on the plane which is halved in the periodic boundary condition.
- We carefully examine the histogram of the internal energy, but do not find any clear signature of bifurcation.