Third boundary of the Shastry-Sutherland Model by Numerical Diagonalization

The Shastry-Sutherland model --- the $S=1/2$ Heisenberg antiferromagnet on the square lattice accompanied by orthogonal dimerized interactions --- is studied by the numerical-diagonalization method. Large-scale calculations provide results for larger clusters that have not been reported yet. The present study successfully captures the phase boundary between the dimer and plaquette-singlet phases and clarifies that the spin gap increases once when the interaction forming the square lattice is increased from the boundary. Our calculations strongly suggest that in addition to the edge of the dimer phase given by $J_{2}/J_{1}\sim 0.675$ and the edge of the N$\acute{\rm e}$el-ordered phase given by $J_{2}/J_{1}\sim 0.76$, there exists a third boundary ratio $J_{2}/J_{1}\sim 0.70$ that divides the intermediate region into two parts, where $J_{1}$ and $J_{2}$ denote dimer and square-lattice interactions, respectively.

It is well known that frustration in magnetic materials enables exotic quantum states to be realized. However, not many such quantum states are obtained in a mathematically rigorous form. The Shastry-Sutherland model 1) is a member of the family of mathematical models in which not all but only some of the eigenstates are exactly obtained. [2][3][4][5][6][7][8][9][10][11] Among such systems, the Shastry-Sutherland model became important after a good candidate material, SrCu 2 (BO 3 ) 2 , was discovered. 12,13) The discovery was followed by extensive theoretical and experimental studies.
Between the region with the exact dimer ground state and the weakly frustrated region with the typical Néel-ordered ground state, the existence of the plaquette-singlet phase was pointed out in Ref. 14. Various approaches [15][16][17][18][19] have theoretically attempted to clarify the behavior of the system in the intermediate region. In the numerical-diagonalization studies, unfortunately, the maximum of the treated system sizes -32 spin sites in Ref. 16 to the best of our knowledge -was not so large. On the other hand, the pressure dependence of the spin gap was observed experimentally. [20][21][22] Among these studies, Ref. 22, employing an electron spin resonance study under high pressure and high field, recently reported the behavior of the spin gap around the phase transition at the edge of the dimer phase. This experimental result becomes a significant motivation to theoretically clarify the behavior of the spin gap around the intermediate region from the numericaldiagonalization calculations for even larger systems.
Under the circumstances, the purpose of the present paper is to report numerical-diagonalization results of even larger systems and to examine the behavior of the system around the intermediate region between the dimer and Néel-ordered regions. The present study provides results for 36-site and 40site systems in addition to smaller systems to deepen our understanding of this system. In particular, numerical diagonalizations of the 40-site system require large-scale parallel calculations in an appropriate supercomputer. We successfully detect the edge of the dimer phase and the edge of the Néelordered phase from the results of the two large sizes. Between the two edges, we additionally detect the third ratio of the boundary dividing the intermediate region into two parts, each state of which has characteristics of a correlation function that are different from each other. The Hamiltonian studied here is given by Here, S i represents the S = 1/2 spin operator at site i. We consider the case of an isotropic interaction in spin space in this study. Site i is assumed to characterize the vertex of the square lattice. The number of spin sites is denoted by N s . The first term of Eq. (1) denotes orthogonal dimer interactions represented by thick solid bonds in Fig. 1. The second term of Eq. (1) represents interactions forming the square lattice represented by thin solid bonds in Fig. 1. We consider that the two interactions between the two spins are antiferromagnetic, namely, J 1 > 0 and J 2 > 0. Energies are measured in units of J 1 ; hereafter, we set J 1 = 1. We denote the ratio J 2 /J 1 by r. Note here that when r = 0, the system is an assembly of isolated dimerized-spin models, whereas the system is reduced to the S = 1/2 Heisenberg antiferromagnet on the ordinary square lattice in the limit r → ∞. We treat finite-size clusters with system size N s under the periodic boundary condition. In this study, N s = 16, 20, 32, 36, and 40 are treated; finite-size clusters are shown in Fig. 1 are tilted. The regular-square clusters help us capture well the two dimensionality of the present system. We carry out our numerical diagonalizations on the basis of the Lanczos algorithm to obtain the lowest energy of H in the subspace belonging to j S z j = M. Note here that the z-axis is taken as the quantized axis of each spin. It is widely believed that numerical-diagonalization calculations are unbiased. Thus, one can obtain reliable information about the system. The energy is denoted by E(N s , M), where M is an integer; in particular, we calculate the cases M = 0 and M = 1 because our attention is focused primarily on the behavior of the spin gap given by Some of the Lanczos diagonalizations were carried out using MPI-parallelized code that was originally developed in the study of Haldane gaps. 23) The usefulness of our program was confirmed in large-scale parallelized calculations. ?, [24][25][26][27][29][30][31][32][33][34] Note here that the largest-scale calculations in this study have been carried out using either the K computer or Oakforest-PACS. Now, let us observe the N s -dependence of ∆/J 1 for some representative cases of r; the results are depicted in Fig. 2. One finds that for r = 1.5, ∆/J 1 significantly decreases as N s is increased. The decreasing behavior of ∆/J 1 is consistent with that in the gapless Néel-ordered phase. Our results for r = 1.5 suggest an almost linear dependence on 1/N s . On the other hand, for r = 0.72 and 0.69, ∆/J 1 decreases with increasing N s for small N s but shows only a very weak N sdependence for large N s . For r = 0.66, ∆/J 1 finally becomes almost constant for all ranges of N s . From these observations, it is considered that clusters with N s = 36 and 40 capture well the behavior of large systems approaching the thermodynamic limit. Therefore, focusing our attention on the results of N s = 36 and 40, we hereafter investigate the behavior of the present system.
Next, let us observe the r-dependence of the spin gap ∆ for finite-size clusters in detail; the results of N s = 36 and 40 are depicted in Fig. 3. First, one finds, in the region up to r ∼ 0.67, that ∆/J 1 gradually decreases as r is increased and that data for the two sizes agree well with each other. The agreement strongly suggests that the system-size dependence has already become weak and that finite-size results almost agree with the corresponding values for the thermodynamic limit. Next, in the region from r ∼ 0.68 to r ∼ 0.70, on the other hand, ∆/J 1 gradually increases with increasing r. The good agreement of data for the two sizes is still maintained, although the r-dependence of whether it increases or decreases has been changed. In the region above r ∼ 0.70, ∆/J 1 decreases once but increases again as r is increased. The upturn of ∆/J 1 is observed for both N s = 36 and 40. The significant characteristic in this region is that there appears a considerable system-size dependence: that is, ∆/J 1 for N s = 40 is smaller than that for N s = 36 at a given r. To find whether the nonzero spin gap exists or is absent in the thermodynamic limit, we will need to carry out further analysis.
To deepen our understanding of the behavior of ∆ showing a complex dependence of decreases and increases, next, let  Fig. 4, therefore, one can understand that the observed changes in the dependence of ∆ are due to the energy-level structure of both E(N s , 0) and E(N s , 1). Next, let us examine the system-size dependence of the spin gap in the region of large r. When r is infinitely large, the system is reduced to the simple square-lattice antiferromagnet, showing that the spin excitation is gapless owing to the existence of the Néel order. As a means of distinguishing whether the system is gapped or gapless, the method of observing the product of the system size and the spin gap is known. This method was successfully used in the study of the plateau --the gap under the magnetic field -at the one-third height of the saturation in the triangular-lattice Heisenberg antiferromagnet with next-nearest-neighbor interactions. 30) The results of this analysis for the present system with N s = 36 and 40 are depicted in Fig. 5. One clearly finds that the results from the two sizes agree with each other in the region down to r ∼ 0.75. When r is further decreased, the results of N = 40 clearly become larger than those of N = 36. The agreement in the behavior of N s ∆/J 1 in the region of large r suggests that the finite-size spin gap in this region exhibits ∆ ∝ 1/N s , which means that the system is gapless. In the region below r ∼ 0.75, on the other hand, the finite-size spin gap does not have the dependence ∆ ∝ 1/N s . Although the N s -dependence of ∆ in the region between r ∼ 0.71 and r ∼ 0.75 is unclear at the present stage, there are two possible scenarios. One is that the system is gapped without any long-range orders. The other is that the system is gapless, but the N s -dependence of ∆ is different from ∆ ∝ 1/N s corresponding to the Néel-ordered phase.
To capture the change from the Néel-ordered phase to the plaquette-singlet phase, let us observe correlation functions in the ground state, namely, S z i S z j ; the results for both N s = 36 and N s = 40 are depicted in Fig. 6. The case r = 1.50 in Fig. 6(b) is a typical one for the Néel-ordered phase. All the results shown by triangles and inversed triangles are positive; this feature is explained by the fact that both i and j for a measured pair are in a common sublattice among the two sublattices of the Néel-ordered state. The results shown by the triangles and the inversed triangles also indicate a gradual decay as the distance is increased. The results shown by the squares indicate alternating signs. This behavior suggests the staggered nature of the Néel-ordered state. Therefore, the characteristics of the Néel-ordered state are well captured in Fig. 6(b). The results in Fig. 6(d), on the other hand, are completely different from those in Fig. 6(b). Among the results shown by the triangles and the inversed triangles, only the shortest-distant datum by the inversed triangle is positive, and the rest are negative. Absolute values of S z i S z j for distances larger than two are very small. These behaviors of the correlation functions are different from those of the Néel-ordered state, but are consistent with those of the plaquette-singlet state. In this state, each plaquette singlet is located at a local square involving the J 1 bond and is the one that has a component of two-spin singlet in diagonal pairs of the square among two possible singlet states of four spins. In the results in Fig. 6(c) for r = 0.73, the pattern of whether S z i S z j is positive or negative is common with Fig. 6(b) and different from Fig. 6(d). A significant difference between Fig. 6(c) and Fig. 6(d) is S z i S z j for the shortest-distant pair shown by the triangle; its r-dependence is depicted in Fig. 6(e) and (f). In Fig. 6(e), the dependence reveals a discontinuity at r ∼ 0.675 for N s = 36 and 40; in Fig. 6(f), another discontinuity appears for N s = 36, which divides the region of r into negative S z i S z j and positive S z i S z j regions. For N s = 40 in Fig. 6(f), S z i S z j changes its sign around r similar to that for the discontinuity of N s = 36, although S z i S z j for N s = 40 is not discontinuous. One possible scenario for the spin state in the larger-ratio region is that, if the system forms plaquette singlets located at a local square involving the J 1 bond, each plaquette singlet is the other singlet state which does not include a component of two-spin singlet in diagonal pairs of the square. Therefore, our calculations suggest that the state for r = 0.73 shows a behavior that is different from that of the state for r = 0.68 and that the behavior changes at common r values for N s = 36 and 40.
To find out whether or not the Néel-type long-range order survives, next, let us observe the r-dependence of correlation functions S z i S z j in detail; the results are depicted in Fig. 7. One finds that as r is decreased down to r ∼ 0.7, S z i S z for the next-nearest-neighbor pair gradually decreases, but its magnitude is not so small. On the other hand, S z i S z j for the pair between the longest distance decreases more rapidly; its magnitude becomes considerably small in the region below r ∼ 0.8. To capture the difference of S z i S z j between the nextnearest-neighbor pair and the longest-distance pair, we examine R cf , defined as the ratio of S z i S z j for N s = 40 divided by the corresponding S z i S z j for N s = 36 presented in Fig. 7; the results are depicted in Fig. 8. One finds that R cf for the longest-distance pair significantly decreases below r ∼ 0.76, whereas the ratio for the next-nearest-neighbor pair is maintained at R cf ∼ 1. This observation suggests that the Néel-type long-range order survives in the region above r ∼ 0.76 and that the order may disappear in the region below r ∼ 0.76, where the Néel-type short-range correlations still survive. As previous estimates for the edge of the Néel-ordered phase, recall r = 0.75 in Ref. 17, and r = 0.765 (15) in Ref. 18; the present result, r ∼ 0.76 for the edge of the region where the Néel-type long-range order definitely exists, agrees well with those previous estimates.
In summary, we have studied the Shastry-Sutherland model by the Lanczos-diagonalization method. The present study has presented diagonalization results for 36-site and 40-site clusters that have not been reported before. Our numerical results have successfully clarified the dependence of the spin gap on the ratio of interactions. Our calculations have successfully captured the edge of the dimer phase to be J 2 /J 1 ∼ 0.675 and the edge of the Néel-ordered phase to be J 2 /J 1 ∼ 0.76. A noteworthy finding is a third specific ratio J 2 /J 1 ∼ 0.70 which divides the intermediate region between the ratios of the two edges into two parts. We have found from observation of correlation functions that the spin state in the smallerratio region and the one in the larger-ratio region are different from each other. The properties of the ground states in the two intermediate regions should be further studied from different viewpoints in future. Such studies would greatly contribute to our fundamental understanding of frustrated magnetism.