# A Study of Highly Frustrated Spin Systems with mixed PEPS in Infinite Honeycomb Lattice

###### Abstract

Highly frustrated spin systems represent a central and challenging problem in condensed mater physics. To this problem, we introduce an algorithm based on mixed projected entangled pair states (m-PEPS), which is a novel type of tensor network. We use the famous Kitaev model on an infinite honeycomb lattice, which can be solved exactly, as a benchmark. With very limited parameters and finite scaling, our calculation results are in good agreement with the exact results, indicating the efficiency of our algorithm. After presenting the benchmark, we investigate the Kitaev-Heisenberg model, which was proposed to describe the effective magnetic momentum interaction in iridate NaIrO which may be used to realize the spin liquid phase. However, our calculations suggest that the gapless spin liquid phase is not robust at the thermodynamic limit, and thus this phenomenon is very difficult to observe.

###### pacs:

02.70.-c, 75.10.Kt, 03.67.Mn, 05.50.+q, 75.10.JmFrustrated spin systems play an essential role in two-dimensional many-body physics and have attracted the attention of condensed matter physicists for decadesbalents (). In particular, some exotic ground states, such as the quantum spin liquid state (which has recently been strongly supported experimentally in the Kagome lattice by Y. Leeexperi () et al), may play a key role in explaining high- superconductivity, as proposed by P. W. Andersonand (). The Kitaev model on a honeycomb latticekitaev (), which can be solved exactly and supports a spin liquid phase, is an important frustrated model in two-dimensional physics. In addition, this proposed model was proposed that it may be realized in a cold atom systemluming (), or in the iridate NaIrO KHmodel (). The exact solution shows that the spin liquid phase is gapless, and a gap opens when a magnetic field, which breaks the time-reversal symmetry, is applied to this model as a perturbation. Unfortunately, general two-dimensional frustrated spin models - even slightly modified Kitaev models such as the Kitaev model perturbed with a magnetic field - cannot be solved exactly.

Due to the lack of analytical methods for general many-body systems, numerical analysis remains the main method for understanding such physics. However, traditional algorithms, such as the quantum Monte Carlo (QMC) method and density matrix renormalization group (DMRG) method, have failed to simulate frustrated systems in two dimensions: QMC suffers from the notorious ‘sign problem’ sign () while DMRG is limited to one dimension and generally lose its power in higher dimensional systemsdmrg1 (); dmrg2 (). Fortunately, the recently developed tensor network (TN) algorithms, such as, the algorithm based on projected entangled pair states (PEPS)PEPS1 (); PEPS2 (); PEPS3 () which is a natural generalization of the DMRG algorithm to higher dimensions, has shown great potential for addressing this problem. The algorithm based on infinite PEPS has been used to obtain a reasonable phase diagram of the frustrated antiferromagnetic Heisenberg model on a checkerboard latticechan (). More recently, a spin liquid phase was claimed near the maximally frustrated region () for the spin antiferromagnetic Heisenberg model on a square lattice for the TN algorithmwangling (); hongcheng (). These results show that the TN algorithms are powerful tools for exploring the frustrated systems and for evaluating some new physics beyond the traditional methods.

Despite the success of the TN algorithm, the study of highly frustrated large systems, such as those represented by the Kitaev model and related models, remains a challenge. The general TN algorithm, such as the algorithm based on PEPS, cannot give the appropriate results for the Kitaev model, because the limited ability of the calculation constrains the bond dimension to a relatively small number (typically ). Thus new methods must be developed to highly frustrated systems with the current calculation abilities. In this letter, we introduce an algorithm based on infinite mPEPS, which is a mixture of the bosonic PEPS (bPEPS)PEPS2 () and the fermionic PEPS (fPEPS)fPEPS1 (); fPEPS2 (). We will demostrate that this novel method can be used to efficiently study highly frustrated systems with a relatively small parameter . We use the Kitaev model as a benchmark and obtain satisfactory results, for both the energy, and the correlations of this model. After the benchmark calculation, we apply this method to the Kitaev-Heisenberg model, which is proposed to describe the effective magnetic momentum interaction model of the iridate NaIrO KHmodel (); comment1 (), to obtain a phase diagram. Our calculations indicate that the gapless spin liquid region is very narrow in the thermodynamic limit and that would thus be very difficult to observe in the iridate NaIrO, even with the Kitaev -Heisenberg model.

Mixed Projected Entangled-Pair States The TN algorithm is a variational method based on TN states with special structures. A general TN state is defined in two steps: first, we describe a configuration of entangled states between virtual particles (we call each entangled state a bond, connected by a line, as shown in Fig.1), which determines the structure and the upper bound of the entanglement of a TN state; second, we define projectors to project the virtual-particle space into real physical space. The different TN algorithms are distinguished by the structure of the TN state, such as PEPS, or the string bond state. The variational parameter space of TN algorithms is determined by the parameters of the projectors, the number of which is dependent on a polynomial of the dimension of the bond , the dimension of the physical space and the number of projectors (as determined by the structure of the lattice and the symmetry of the state). For the special TN state, PEPS, the configuration of the bonds is completely determined by the lattice where each edge of the lattice corresponds to a bond, and each site has an on-site projector. Without a loss of generality, we use the honeycomb lattice as an example. In general, the entangled pair can be written as: , where is the state of virtual particle with dimension and denotes an edge connecting sites and in the lattice. Algebraically, the state can be represented by local creation operators from the vacuum through the second quantization method. In the original PEPS case, the entangled pair state can be represented as bosonic operator as: .

Similarly, the projector of each site on a honeycomb lattice can also be represented as bosonic operators as: , where is the coordinate of the site, is the bosonic creation operator of the physical particle, and ,, are the bosonic annihilation operators of the virtual particles connected to edges , and , respectively (see Fig.1). In this case is the vacuum of the physical space, and is the vacuum of the virtual space. Using the representation of the entangled pairs and projectors, the PEPS state can be represented as: , where all of the virtual particles are projected as physical particles.

In general, the virtual particles in an EPR pair are not necessarily bosons. In solving the many-body physics of fermionic systems, it is natural to extend the bPEPS formulism to the fPEPS formulism. This can be completed by replacing the bosonic operators in the entangled pairs and projectors, as shown in Eq.(1) with fermionic operators. In the fPEPS formulism, the parity of the projectors at different sites of the lattice, which determine the exchange character of the projectors, should be the same even. This requirement guarantees that the fPEPS is well defined, that is, the order of the projectors is not important for the state beyond the global phase.

Here we focus on a frustrated spin system, which exhibits bosonic statistics in physical space. Thus, we propose an algorithm based on the mixed PEPS (mPEPS), which introduces fermionic statistics to the virtual-particle space (same as fPEPS) and bosonic statistics in the physical-particle space (different from fPEPS). Therefore, the projectors of the mPEPS should be expressed in the second quantized form as: , where satisfies: and ,, satisfy: . To validate this definition, the parity of the fermionic operators in the projector should also be same for every site (for example, it could be even). Although there is some overlap, the state spaces of the PEPS, fPEPS and the mPEPS are expected to be different, even for the same parameters.

After the defining mPEPS, the ground state of a system in variational space can be found by the imaginary evolution method which is performed in the same manner as in the PEPS method; in particularly, we use the simple update method to approximate the environment in the same manner as PEPS in two dimensionsxiangtao (). When the ground state is stable under the imaginary evolution, the physical quantities of the system can be calculated through a complex contraction process (for more details, see the Supplementary Information). In the contraction, the bosonic and fermionic operators must be carefully handled(see in the Supplementary Information).

Kitaev Model as the Benchmark Here, we use the Kitaev model, which is a highly frustrated spin model in two dimensions, as a benchmark of the mPEPS algorithm. The Kitaev model is defined on an infinite honeycomb latticekitaev (), and its interaction along different directions (see Fig.1) vary strongly as: . This model can be exactly solved with two resultant phases: the gapped phase and the gapless spin liquid. There are two types of quasi-particles in the gapped phase: fermions and vortices. The vortices associated with a gauge field and a factor will appear when the fermion moves around the vortex. A gap opens when a magnetic field is applied to the gapless spin liquid phase. The non-abelian anyon, which may play an important role in topological quantum computation, will also appear in this phase. In our calculation, we normally set the parameters of the model to satisfy and , which gives a line connecting a vertex to the middle point of the corresponding edge (see the red dashed line in the inset of Fig.2). After conducting the stability tests of the algorithm with respect to the parameters (particularly, the truncation parameter in the contraction process, see details in the Supplementary Information), we compare the energy of the system along the chosen line with the exact solution. Result show good agreement (see the Supplementary Information). In addition, the nearest neighbor(NN) correlation which corresponds to along the line is also shown in the Fig.2. The black line shows the exact value of this correlation as derived inbaskaran (), which is in good agreement with our calculation (red dots). The inset shows the , which exhibits a sharp transition point near 0.46 and indicates a phase transition in this model, although it slightly deviates from the exact result of 0.5.

We must emphasize that the virtual dimension of the system is rather limited by our calculation ability. We can only calculate results for the dimension up to 8. And we conduct the finite scaling for this limited dimension in this model, and the points shown in Fig.2 were obtained with those finite scaling. To compare our findings with the traditional PEPS case, we also present the results from the PEPS with the same parameters and finite scaling employed in Fig.2 (blue triangles). It is clear that the present method improves the results obtain with limited calculation ability. More importantly, the PEPS method does not give an indication of the phase transition, while our mPEPS technique can. However, in simpler cases, such as two-dimensional Ising model with a transverse magnetic field on the honeycomb lattice, the mPEPS and PEPS methods give nearly the equivalent results.

Analytical resultkitaev () suggest that the gapless spin liquid will exhibit a gap upon the addition of a magnetic field in the direction ( direction) as a perturbation, which breaks the time-reversal symmetry. Unfortunately, the Kitaev model with magnetic field cannot be solved exactly. We can only study this transition by numerical methods. Our mPEPS calculations show this transition directly by where is the magnitude of the magnetic field (see Fig.3), and these calculations strongly suggest a first order phase transition occurring at very low magnetic fields, .

Application to the Kitaev-Heisenberg Model Having verified the mPEPS method in the Kitaev model, we extend the application of this method to the Kitaev-Heisenberg model. The Hamiltonian of this model can be written as: , where is the standard Heisenberg interaction: and the is Kitaev interaction. With the Khaliullin transformationKHmodel (), which performs rotations on different sites, this model can be exactly solved at the point and . When , the ground state of the transformed system can be exactly solved, and is found to be a ferromagnetic state. With the inverse Khaliullin transformation, the ferromagnetic state will be transformed to the stripy anti-ferromagnetic state. When , this model is equivalent to the Kitaev model, with , whose ground state can be exactly solved as a gapless spin liquid state. For the parameters beyond and , the model can not be exactly solved and can only be examined by numerical calculations. Former studies have employed several numerical methods, such as exact diagonalizationKHmodel (); Schaffer () and DMRGJiang (), and have found three different phases: the Neel order phase, the stripy anti-ferromagnetic phase and the spin liquid phase. These previous calculations, which were based on the systems of relatively small finite sizeKHmodel () or quasi one-dimensionalJiang (), showed that phase transition points occur near and at zero temperature.

We applied the mPEPS method to this model on an infinite two-dimensional lattice which naturally satisfies the thermodynamic limit. We directly calculate the spin configuration of the ground state with which is indeed a ferromagnetic state, as predicted by the exact results. As with the benchmark calculation of the Kitaev model with (corresponding to ), here we also give , which can be used to determine the phase transition point, as shown in Fig.4.

In our infinite system with six-site periodicity, the second transition point is dramatically different from the former result(), while the first transition point, at , is slightly different from . To confirm our calculation, we also determine the NN correlation , which is used as an indication of the phase transition in referenceKHmodel (). The results also support the occurrence of a phase transition at the same position. In our calculations, the second phase transition point occurs at for , and at for . This findings suggest that the transition point will tend to with the finite scaling for . Our results indicate that the gapless spin liquid phase will be destroyed by the Heisenberg perturbation very quickly, and thus this phase is not robust for this type (Heisenberg) of perturbation. The direct implication of our results indicates that it would be very difficult to observe a Kitaev-type gapless spin liquid in iridate NaIrO. The findings of a decreasing spin liquid phase in the infinite system (at the thermodynamics limit), compared to the finite exact diagonalization and DMRG in quasi one-dimension with cylinder boundary condition, is not so surprising. The size of the system plays a subtle role in the spin liquid phase. Z. Meng et al have reportedZYMeng () a spin liquid phase between the semimetal phase and an antiferromagnetic Mott insulator in the Fermi-Hubbard model on the honeycomb lattice, with the inclusion of 648 sites, by using projective QMC method. However, when S. Sorella et aljapanese () reconsidered this problem by including up to 2592 sites, the evidence of spin liquid was lost. In our case, the infinite lattice may play a similar role.

To summarize, we have introduced a novel type of TN (mPEPS) algorithm to attack the highly frustrated spin systems. Our method is successful in simulating the Kitaev model, with respect to both the energy and the correlations of the ground state. Having established this method, we applied it to the Kitaev-Heisenberg model. Our calculations show that the spin liquid region is very narrow, and thus this state would be very difficult to observe experimentally. The success of our method opens the possibility of considering virtual particles in a tensor network with more complex statistical characteristics, such as the fractional statistics of anyons. These new statistics may be convenient for studying the system which excitations have the similar statistics.

We acknowledge support from the Chinese National Fundamental Research Program 2011CB921200 and the Fundamental Research Funds for Central Universities Nos. WK2470000006, WK2470000004, WJ2470000007 and NSFC11105135.

## References

- (1) P.W. Anderson, Science 235, 1196 (1987).
- (2) Leon Balents, Nature 464,(2010).
- (3) Tian-Heng Han et al Nature (2012).
- (4) A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
- (5) L.-M. Duan, E. Demler, M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
- (6) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- (7) U. Schollwock, Rev. Mod. Phys. 77, 259 (2005).
- (8) C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81, 052338 (2010).
- (9) P. Corboz, R. Orus, B. Bauer, G. Vidal, Phys. Rev. B 81, 165104 (2010).
- (10) Actually, the Kitaev-Heisenberg model is not the only model proposed for describing the effective magnetic momentum interaction model of iridate NaIrO. Other models have proposed. Such as that described in: Y. Singh et al, Phys. Rev. Lett. 108, 127203 (2012). In this work we only consider the Kitaev-Heisenberg model.
- (11) G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev. Lett. 98, 247201 (2007)
- (12) L. Wang, D. Poilblanc, Z.-C. Gu, X.-G. Wen, and F. Verstraete, arXiv:1301.4492.
- (13) H. C. Jiang, Hong Yao and L. Balents, Phys. Rev. B 86,024424(2012).
- (14) Y. H. Chan, Y.J. Han, L. M. Duan, Phys. Rev. B 84, 224407(2012).
- (15) E. Y. Loh, et al Phys. Rev. B 41, 9301(1990)
- (16) J. Jordan, R. Orus, G. Vidal, F. Verstraete and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
- (17) F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066.
- (18) F. Verstraete, J.I. Cirac, V. Murg, Adv. Phys. 57,143 (2008).
- (19) J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 105, 027204 (2010)
- (20) H.-C. Jiang, Z.-C. Gu, X.-L. Qi, and S. Trebst, Phys. Rev. B 83, 245104 (2011)
- (21) R. Schaffer, S. Bhattacharjee and Y. B. Kim, arXiv:1206.5814.
- (22) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad and A. Muramatsu, Nature 464, 847(2010)
- (23) S. Sorella, Y. Otsuka and S. Yunoki, Sci. Reports 2, 992 (2012).
- (24) H. C. Jiang, Z. Y. Weng, and T. Xiang, PRL 101, 090603 (2008).

## .1 Supplementary

### .1.1 Imaginary Time Evolution

After chosen the variational space of the state, i.e., mPEPS with given parameter and , we use the imaginary time evolution methodPEPS1 (); PEPS2 (); PEPS3 () to find the ground state. The imaginary time evolution could be summarized as the following formula:

The is the initial state, which can be chosen randomly, after long enough imaginary time evolution, the final normalized state will be close to the ground state with well controlled precision, as long as the initial state is not orthogonal with the ground state (with zero measurement to chosen these states). Generally, we can only deal with the local evolution operators (the time evolution of a system is a global operator), and the Trotter expansion can replace the global operators by local operators with controllable precisionPEPS1 (); PEPS2 (); PEPS3 (). For example, in Kitaev model, Hamiltonian could be split as , where is the part of Hamiltonian acting on -bond. With the help of Trotter expansion, we could split the evolution operator in second order as:

With the Trotter expansion, it seems that we just need to act the local evolution operators on the initial state with the given order and obtain the ground state finally. However, the dimension of the bond will increase exponentially with the evolution time , if local operators are applied without approximation. Therefore, as in the general TN algorithm, approximations should be made after each step of evolution. One of the simplest way to cut the tensors back to the initial parameter space is the singular value decomposition (SVD), which maintain the important part of the tensor (see Fig.A1). It is straightforward to use SVD in the case of PEPS, while it needs extra operations in the case of mPEPS (also in fPEPS) with fermionic statistics. Since the tensors in mPEPS are set with even parity, the matrix in the SVD process can be re-indexed as block diagonal matrix according to parity of the tensor’s indices (one block with odd parity and the other with even parity), then we can apply the SVD and cut off in each of the block.

### .1.2 Contract

After the imaginary time evolution, we suppose that we have found the ground state. Then the key task is to find the values of some physical quantities which can be compared with the experiment. This task is realized by the contract processes, which including the contract of tensors and entangled pairs. The physical quantities in a given ground state can be expressed as

For simplicity, we use as an example to show how to contract to calculate the physical quantities. The calculation of physical quantities is similar. Suppose we have a 2D honeycomb lattice and each site is labeled with two indices i,j. We have a state in the form of projected entangled-pair State:

where and denotes the edge e, linking site and site . Normally, translational invariance of tensors will be imposed in infinite lattice. Actually, we have the corresponding bra vector as

where . The product order is reversed which will make a difference in the contract procedure below. The scalar product of bra vector and ket vector is:

Because we are considering mPEPS where tensors have even parity, it is free to change the orders. We would like to rearrange the order as:

The product of will contract over physical indices and we symbolize it as .

where repeated indices mean contraction. Notice that is not just the contraction of and , but with extra minus signs resulting from commutation of operators.

Therefore, the scalar product is

What we should do next is to contract over all virtual indices, ie. to product all the tensors by EPRs. Even though we consider here is a infinite lattice, we suppose there are boundaries in very remote areas. We contract row by row many times from boundaries until the boundary tensors are stable. After we contract the boundary with a row of tensors, the boundary tensors will get thicker, like the case of imaginary time evolution, and we need to do approximations to cut them back to their initial shape. Normally we use a specific form of boundary tensors to make the approximation, for the sake of simplicity (see Fig.A2 for details).

The main tool used to make this approximation is still SVD: We fix two tensors, and trace over left and right parts of boundary tensors to gain the left and right environments of the two fixed tensors. Now we contract over left environment tensor, fixed tensors, right environment tensors together to become a bigger tensor. We use SVD to cut the bigger tensor, then cancel the environment effect by contracting inverse tensors of the left and right environments. It is extremely obvious if you write down the approximation judgement explicitly:

where is the old boundary contracting with a row of tensors, B is the new boundary with only two tensors different from the old boundary. G is the unit cell of . are perspectively the left and right part of G in the row of . L is just and R is just . The above equation is just binary form of and . Obviously, reaches its minimum when best approximates . See also Fig.A3.

This example of showed our method for contract. It is similar to calculate average value of local operators , except that the sites where operators lie should be kept until final contract.

There is one more point worth mentioning: When we use mPEPS, we choose the virtual annihilation and creation operators to be anti-commutative. So we must be extremely careful about the minus sign when we contract over virtual indices. Similarly and more easily to overlook, in SVD, minus signs also appear when we change the order of virtual indices in order to obtain the original shape of tensors. It could be explained as equations below. We take a simplified version as an example. T is a big tensor which needs to be split by SVD.

The first equation is normal SVD; the second is to absorb singular values; and rest equations are to extract EPR from tensors. EPR is , correspondingly , and so as their Hermitian conjugate. In the equations, you can see that a minus sign show up. Note that in the fifth equation is replaced by since it is not zero only when and . p(*) is the function that maps quantum numbers to parity. U and V is the tensors that and absorb square of the singular values respectively. is the tensor after U absorbs .

### .1.3 Additional Calculation of the Kitaev Model

In order to shows the stability of our results, we increase the dimension of the truncation n the contract from to . We can see in the following figure that when with , the NN correlation of the system is relatively stable.

With the stable algorithm, we calculate the energy of the Kitaev model which is well agree with the exact energy

Since the limit of our calculation, we should do some finite scaling. However, the finite scaling of D is not clear and well defined. Here we use linear of polynomial fit with , as showed below.

### .1.4 Khaliullin Transformation for Kitaev-Heisenberg Model

This transformation first proposed by Khaliullin and Chaloupka not only makes the anti-ferromagnetic state to ferromagnetic state, but also makes Kitaev-Heisenberg model exactly solvable at . This transformation requires different spins in the lattice to rotate about different axis. To be specific, we first choose a set of spins which are positioned on third nearest neighbor sites at opposite corners of the hexagons throughout the lattice, and hold these spins fixed. We next rotate the three nearest spins by about the spin axis corresponding to the bond which connects it to the fixed spin. The effect of this operation therefore transforms the Heisenberg Hamiltonian to a mixture of Heisenberg Hamiltonian and Kitaev Model, and Kitaev Hamiltonian invariant.

Thus the Kitaev-Heisenberg model is transformed as:

(A1) |

It is obvious to find out in the transformed Hamiltonian, when , it is just Heisenberg Hamiltonian with a negative coupling coefficient. So the state at after the transformation is a ferromagnetic state. We can operate an inverse transformation to send the ferromagnetic state back to an anti-ferromagnetic state. As we have explained before, the anti-ferromagnetic state is harder to calculate, since stable results are not easily required. So in the letter, we just use the Hamiltonian after the transformation to discuss the quantum phase diagram of Kitaev-Heisenberg model.

In order to weaken the influence by Heisenberg interaction and make the spin liquid region broader and easier to see, we re-parameterize the Kitaev-Heisenberg model to:

(A2) |

The new parameter and the old one can be related by a simple formula, if is proportional to to the same extent. Besides, the new Hamiltonian Eq.A2 preserves the same behaviors as Eq.A1 at and . In the paper, we discuss Eq.A2, and note that the spin liquid region of Eq.A1 is actually much much narrower. The region of spin liquids is , as the part below suggests. Correspondingly, the real parameter region of spin liquids before re-parameterization are around . In the paper, we just use the new Hamiltonian Eq.A2 to further our discussion.

The configuration of the ground state of the Kitaev-Heisenberg model with can be calculated by our PEPS method as Fig.A8