Signatures of Quantum Spin Liquid in Kitaev-like Frustrated Magnets

Signatures of Quantum Spin Liquid in Kitaev-like Frustrated Magnets

Abstract

Motivated by recent experiments on -RuCl, we investigate a possible quantum spin liquid ground state of the honeycomb-lattice spin model with bond-dependent interactions. We consider the model, where and represent the Kitaev and symmetric-anisotropic interactions between spin-1/2 moments on the honeycomb lattice. Using the density matrix renormalization group (DMRG), we provide compelling evidence for the existence of quantum spin liquid phases in an extended region of the phase diagram. In particular, we use transfer matrix spectra to show the evolution of two-particle excitations with well-defined two-dimensional dispersion, which is a strong signature of quantum spin liquid. These results are compared with predictions from Majorana mean-field theory and used to infer the quasiparticle excitation spectra. Further, we compute the dynamical structure factor using finite size cluster computations and show that the results resemble the scattering continuum seen in neutron scattering experiments on -RuCl. We discuss these results in light of recent and future experiments.

1Introduction

One of the hallmarks of quantum spin liquid is the existence of fractionalized excitations[1]. While it is generally difficult to detect the dispersion of single-particle excitations in fractionalized systems, the information about two-particle and other multi-particle excitations is contained in the dynamical spin structure factor measured in inelastic neutron scattering. If the ground state is a quantum spin liquid, the lower boundary of multi-particle continuum should have a well-defined dispersion. Many candidate materials for quantum spin liquid, however, do not allow such scattering experiment due to unavailability of large single crystal. Possible experiments are, therefore, quite often limited to thermodynamic measurements such as specific heat and susceptibility, as well as thermal transport. Hence it is difficult to identify smoking-gun evidence for quantum spin liquid.

In this context, a recent inelastic neutron scattering experiment on -RuCl provides valuable information about multi-particle continuum in a putative quantum spin liquid material[2]. -RuCl is one of the candidate materials that support the Kitaev interaction[4] between pseudo-spin moments, which is a bond-dependent frustrated Ising interaction and arises from the combination of correlation effects and spin-orbit coupling[5] in Ru ions (see Refs. for recent reviews). When only the Kitaev interaction is present, such a model can be exactly solved[4] and the ground state is a quantum spin liquid. On the other hand, other interactions are generally present and they may drive a transition to a magnetically ordered state[8].

The compound -RuCl orders magnetically at low temperatures[12], possibly due to the existence of other interactions mentioned above. In spite of this, the dynamical structure factor shows a continuum of excitations at high energies both below and above the ordering temperature, which may be related to a nearby quantum spin liquid[15]. Elaborate ab initio computations indicate that the dominant exchange interactions are the Kitaev and symmetric-anisotropic interactions with small third neighbor Heisenberg interaction[9]. A previous theoretical work[16] of exact diagonalization (ED) on finite size clusters suggest that the model may host quantum spin liquid phases in an extended region of the phase diagram and small would drive a transition to the zig-zag order that is seen in the experiment. Such a study is naturally subject to finite size effect and it is still difficult to nail down the precise nature of the ground state. In fact, the ED study mentioned above introduce a spatial anisotropy in the Kitaev interaction to avoid strong finite size effect.

In this work, we investigate the model using then infinite density matrix renormalization group[17](iDMRG), where the system is placed on an infinite cylinder. We first study the ground state energy, entanglement entropy, magnetization, and static structure factor for the model. Similar to the previous ED study, we find quantum paramagnetic ground states for ferro-like Kitaev interaction for an extended region in the phase diagram. It is also shown that the entanglement entropy in this region remains relatively high. In comparison, the entanglement entropy of a magnetically ordered state for anti-ferro-like Kitaev interaction is quite small.

In order to test the existence of coherent excitations, we compute the eigenvalues of the transfer matrix (TM) in the matrix product wavefunction. Using the mapping[20] between the complex eigenvalues and low-energy excitation spectra as a function of momentum, we identify the lower edge of the multi-particle excitation continuum as a function of momentum[21]. From this, we find that the lower edge of the excitation spectra, which also corresponds to the longest correlation length in the system, moves coherently in momentum space in a well-defined fashion as a function of . This alone tells us that these states are non-trivial and likely corresponds to quantum spin liquid phases. Assuming that the lowest momentum-dependent eigenvalues correspond to the lower boundaries of the two-particle continuum in quantum spin liquid, we compare the results of the transfer matrix computations and the two-particle spectra obtained in the Majorana mean-field theory for the model. It is found that the transfer matrix results can be interpreted as coming from the convolution of two single-particle spectra in the Majorana mean-field theory for some representative values of . This suggests that the mean-field picture can capture some essential features of the two-particle excitation spectrum seen in the transfer matrix eigenvalues.

As pointed out in earlier works[22], the dynamical spin structure factor in the Kitaev model involves not just two-particle spectra, but also the flux degrees of freedom. In fact, this makes it difficult to extract the information about the single-particle dispersion from the dynamical spin structure factor as the flux degrees of freedom may still play important roles even in the presence of the additional interaction. In contrast, the transfer matrix eigenvalues computed here are directly related to the convolution of the single-particle spectra.

In order to make a more direct contact with scattering experiment, we compute the dynamical structure factor on a 24-site cluster. It is found that the low energy dynamical structure factor exhibit a star-like feature in momentum space, just like what is seen in recent neutron scattering experiments[3] on -RuCl. We also demonstrate that these results are consistent with the equal-time spin structure factor computed in the Majorana mean-field theory. We have further confirmed that the magnetic-field dependence of the magnetization in the 24-site cluster is consistent with the experimental results[25] on -RuCl. Taken together, all of the results above suggest that -RuCl may be close to a quantum spin liquid described by the model with ferro-like Kitaev interaction.

The remainder of the paper is organized as follows. We present the model and iDMRG results in Section 2. In Section 3 we present and discuss the transfer matrix (TM) spectrum. A Majorana based MFT is presented in Section 4, while results of the ED calculation of the dynamic structure factor appear in section Section 5. Details of the iDMRG and ED calculations, as well as additional results are given in the appendices.

2DMRG study of the model

Figure 1: Ground state energy density, E_{GS}, of the K-\Gamma model, and the corresponding entanglement entropy, S_E, for a bipartition of the cylinder into two infinite cylinders, as determined using iDMRG for cylinders with circumference L=6. Kinks in E_{GS} are indications of first order phase transitions. Inset shows the cylinder’s three-unit-cell circumference.
Figure 1: Ground state energy density, , of the model, and the corresponding entanglement entropy, , for a bipartition of the cylinder into two infinite cylinders, as determined using iDMRG for cylinders with circumference . Kinks in are indications of first order phase transitions. Inset shows the cylinder’s three-unit-cell circumference.
Figure 2: Staggered magnetization (top) and static spin structure factor (bottom), calculated using iDMRG. Inset: Brillouin zone with labeled positions of symmetry points.
Figure 2: Staggered magnetization (top) and static spin structure factor (bottom), calculated using iDMRG. Inset: Brillouin zone with labeled positions of symmetry points.
Figure 3: Kitaev limit transfer matrix spectrum for L=6. Top: isotropic Kitaev model - K_x = K_y = K_z = -1. Bottom: anisotropic Kitaev model - K_x = -1.1135, K_y = 1, K_z = -1.2, such that the gapless nodes remain on the cylinder’s allowed momentum cuts.
Figure 3: Kitaev limit transfer matrix spectrum for . Top: isotropic Kitaev model - . Bottom: anisotropic Kitaev model - , such that the gapless nodes remain on the cylinder’s allowed momentum cuts.
Figure 4: (a) Single Majorana fermion spectrum \varepsilon_c({\bf k}) for the isotropic Kitaev Kitaev model, K_x=K_y=K_z=-1, on a cylinder with a three unit cell circumference. (b) Corresponding minimum energy for two-particle excitations, \Omega_{\rm
      min.}({\bf q}). (c) Majorana spectrum, and (d) corresponding two particle excitation edge, for anisotropic Kitaev model, K_x =
    -1.1135, K_y = 1, K_z = -1.2.
Figure 4: (a) Single Majorana fermion spectrum for the isotropic Kitaev Kitaev model, , on a cylinder with a three unit cell circumference. (b) Corresponding minimum energy for two-particle excitations, . (c) Majorana spectrum, and (d) corresponding two particle excitation edge, for anisotropic Kitaev model, .
Figure 5: Transfer matrix spectrum for different \phi. Plotted is the spectrum E(k) = - \log |\lambda| with each point corresponding to a single eigenvalue \lambda, where \lambda =
    |\lambda| e^{i\eta} and \eta is identified with the momentum k_x along the cylinder. k_y denotes the transverse momentum obtained as a quantum number with respect to translation along the cylinder. See App. A for more details.
Figure 5: Transfer matrix spectrum for different . Plotted is the spectrum with each point corresponding to a single eigenvalue , where and is identified with the momentum along the cylinder. denotes the transverse momentum obtained as a quantum number with respect to translation along the cylinder. See App. A for more details.

The Hamiltonian is given by

where are spin- operators at site of a honeycomb lattice. Unless otherwise noted, we consider here the isotropic case, . Throughout the following, and are parameterized using such that and . We use the iDMRG method with bond dimensions of up to to obtain the ground state of this model on a narrow infinite cylinder with a three unit cell circumference (). In Figure 1 we plot the ground state energy density, , as a function of the parameter . Starting from the ferromagnetic Kitaev limit, , evolves smoothly through the limit, . A discontinuity appears at and again slightly before the anti-ferromagnetic Kitaev limit, . The two discontinuities are associated with a transition into, and out of, a magnetically ordered vortex state, which becomes an exact product state for . This is evident in a plot of the entanglement entropy, , also in Fig. Figure 1, showing a vanishing at this point. Notice that the entanglement entropy remains as high as that of the ferro-like Kitaev limit in the entire region of . Furthermore iDMRG shows a finite magnetization in , as well as an enhanced spin structure factor, Figure 2, all consistent with a magnetically ordered phase. The main question we want to address here is: what is the nature of the ground state outside of the magnetically ordered state? The large entanglement and lack of magnetic order suggest that the ground state in this region is occupied by quantum spin liquid phases. However, a small discontinuity at in both the entanglement entropy and the spin structure factor, raises the question whether there exists a subtle transition between different kinds of spin liquid phases. To address this issue, and to gain insight into the low energy physics of the model as is tuned from the Kitaev to limits, we turn to a detailed examination of the transfer matrix spectrum, obtained from the ground state MPS.

3Transfer matrix spectrum

3.1The Kitaev limit

We begin by analyzing the transfer matrix spectrum, , of the pure Kitaev model, shown in Figure 3, as a function of the momentum along the cylinder, . We were also able to resolve the transverse momentum , which are depicted in the figure by different colors (see Figure 7). Here we use the demonstrated correspondence[20] between the complex eigenvalues of the transfer matrix and the excitation spectrum, . Namely, given a TM eigenvalue , the corresponding momentum (along the infinite dimension) is given by , while the corresponding energy is given by (see Appendix A for details). The Kitaev model is exactly solvable in terms of Majorana fermions, and therefore it is possible to readily identify the features in Figure 3 with the known Majorana excitations. The most prominent feature of the Majorana spectrum, , in the Kitaev model is the existence of two gapless Dirac nodes at the corners of the Brillouin zone, and (see Figure 4a). Thus, single Majorana excitations can account for one of the gapless points in Figure 3a (green pillar), but not for the other. To account for the gapless point at , we must assume that, at least, two-particle excitations, , show up in the TM spectrum as well (Figure 4b). Note, however, that in the isotropic Kitaev model, the locations of the gapless single particle excitations in momentum space coincide with some locations of the gapless two-particle excitations. To verify that the TM spectrum can also depict single particle excitations, we studied the TM for an anisotropic Kitaev model, where the gapless Dirac nodes move from their original positions at and , but with an anisotropy such that the nodes still sit on the cylinder’s allowed momenta cuts. Figure 3b shows that the gapless single particle and two particle excitations at split, appearing exactly where they are expected according to the Kitaev solution, Figure 4c,d, clearly show that both show up in the TM spectrum, as stated in Ref. . As a minimal interpretation of the TM spectrum we therefore suggest that the pillars indicate the positions of minima in the single and two-Majorana spectra. Thus, returning to the isotropic case, Figure 3a, and in accordance with the Kitaev solution, the pillar at is associated with a gapped minimum in the two-particle spectrum, whereas the pillar at is associated with a gapped minimum in the single particle spectrum.

3.2The model

Figure 6: Energy density per bond, as obtained using iDMRG, for systems with a three (L=6) and six (L=12) unit cell circumference.
Figure 6: Energy density per bond, as obtained using iDMRG, for systems with a three () and six () unit cell circumference.

Moving away from the exactly solvable Kitaev limit, we now turn to analyze the TM spectrum of the model, Figure 5, which shows the transfer matrix spectrum, , for various values of . Similarly to the Kitaev limit, minima in the continuum of excitations are clearly identified at , , , and . All, however, are gapped. This can be understood in the context of Majorana fermions by noting that the cylindrical geometry breaks the symmetry between bonds and bonds, which in turn can lead, for , to anisotropic hopping amplitudes, and the gapping out of the fermions. To corroborate this point, Figure 6 depicts the energy density per bond as a function of , displaying that indeed the symmetry between bonds is broken for .

Several additional minima continuously move as is tuned from the Kitaev limit, , to the limit at . For example, close to the Kitaev limit, , there is a minimum near , but as increases, it moves towards at , and then turns back, crossing near . Similar evolutions can be seen for minima at . It is natural to suggest that these additional minima are associated with single particle excitations, since under some conditions they seem to have the lowest energy. However, their positions seem to be related to each-other by addition of or . For example, for , there is a minimum at , which by pairing with a soft excitation at , and inversion, , gives a two-particle minimum at , which is also observed in the TM spectrum. The minimum at is similarly obtained by pairing with a soft excitation at followed by inversion. Therefore, the positions of these minima do not rule out the possibility that for finite , only two-particle excitations are observed in the TM spectrum. Figure 7 shows the coinciding positions of the soft single and two-particle excitations, for and .

Whether we interpret the TM spectrum as being associated only with two-particle excitations, or with the single-particle spectrum as well, we can still reach the following conclusions: (i) The TM excitation spectrum evolves without any qualitative change through , suggesting that the fundamental nature of the ground state is likely to remain the same. Supporting this conclusion is the fact that at there is a sudden increase in bond anisotropy (see Figure 6), which is a result of the cylindrical geometry with finite circumference. Hence the discontinuity can be interpreted as a “Fermi surface topology” change in quasiparticle spectrum. (ii) The positions of the excitation continuum minima, and their evolution as a function of , strongly indicate that the paramagnetic phase of the model harbours coherent excitations commonly associated with quantum spin liquids. Thus, it is reasonable to conclude that in the entire region , the model has a quantum spin liquid ground state, which shares basic properties with the ground state of the ferromagnetic Kitaev model.

Figure 7: Allowed momentum cuts for the cylindrical geometry. Symbols indicate the coinciding positions of the soft single and two-particle excitations, as deduced from the TM spectrum. \bullet indicate the leading soft modes at K,K' and M points. Additional soft modes are plotted for \phi=0.03\pi\,(\blacktriangle) and \phi=0.2\pi\,(\blacksquare). As \phi is tuned between these values, we expect that these minima smoothly move from the squares (\blacksquare) to the triangles (\blacktriangle).
Figure 7: Allowed momentum cuts for the cylindrical geometry. Symbols indicate the coinciding positions of the soft single and two-particle excitations, as deduced from the TM spectrum. indicate the leading soft modes at and points. Additional soft modes are plotted for and . As is tuned between these values, we expect that these minima smoothly move from the squares () to the triangles ().

4Majorana mean-field theory

Motivated by the iDMRG results of the previous section, we would like to formulate a Fermionic mean-field theory which closely resembles the exact solution of the Kitaev model. Therefore, following Kitaev[4], we replace the spin operators in the Hamiltonian with products of Majorana fermion operators, ,

where for a -type bond,

Similar definitions follow for and -type bonds. Here, the Majorana fermion operators are normalized such that and . The physical Hilbert space of the spin Hamiltonian is then obtained by projecting the Majorana Hamiltonian onto the subspace of states which obey . Within a mean-field approach, we can approximate with

where the fields and obey the mean-field self-consistency equations on each bond,

Given the ground state of , it is possible to construct an approximate ground state for by projection onto the physical Hilbert space, .

Figure 8: Band structure of c (red) and b (blue) Majorana fermions, plotted along high symmetry lines of the Brillouin zone, for several values of \phi in the two dimensional thermodynamic limit. The flat bands for \phi=0 are three-fold degenerate; for \phi/\pi=0.5, the lower energy flat bands are two-fold degenerate. When \phi/\pi\sim 0.2, there is a finite density of zero energy b fermion states, around the \rm K (and \Gamma) points in the Brillouin zone. At \phi/\pi=0.4 the b fermion bands are still dispersive, but gapped.
Figure 8: Band structure of (red) and (blue) Majorana fermions, plotted along high symmetry lines of the Brillouin zone, for several values of in the two dimensional thermodynamic limit. The flat bands for are three-fold degenerate; for , the lower energy flat bands are two-fold degenerate. When , there is a finite density of zero energy fermion states, around the (and ) points in the Brillouin zone. At the fermion bands are still dispersive, but gapped.

It is straightforward to obtain a uniform, Z-flux-free, solution of Equations (Equation 5) and ( ?), in the two-dimensional thermodynamic limit. In the following we use the convention that in etc., the subscript indicates a site on the odd sublattice and a site on the even sublattice. Assuming that on all bonds, we obtain , which is independent of . Similarly, is independent of , but it does depend on the ratio . The mean-field ground state energy per bond is given by . By solving , it is possible to obtain the Majorana fermion spectrum, shown in Figure 8 along high symmetry lines in the Brillouin zone, for several values of . In the Kitaev limit, , one finds a single dispersing -fermion band, with Dirac nodes at the K points of the Brillouin zone, and whose band width is set by . Three flat bands describe the -fermions, which are localized on the bonds. For , one obtains a similarly dispersing -fermion band, with band width , and three flat bands which describe the -fermions localized on the hexagons. When both and are nonzero, the -fermions are dispersive, and become gapless in the parameter range .

Figure 9: Left: b (solid) and c (dashed) Majorana fermion mean-field spectrum plotted along the three momentum cuts allowed on a cylinder with a three unit cell circumference (see fig. ). Right: Corresponding two-Majorana b+c (solid) and c+c (dashed) spectrum.
Figure 9: Left: (solid) and (dashed) Majorana fermion mean-field spectrum plotted along the three momentum cuts allowed on a cylinder with a three unit cell circumference (see fig. ). Right: Corresponding two-Majorana (solid) and (dashed) spectrum.

We have suggested that the TM spectrum can be associated with single and two-particle continua of fractionalized excitations. The Majorana MFT can be conveniently used to test this idea. Thus, to compare the two approaches, we now apply the MFT to the cylindrical geometry, considered above using iDMRG. As in the iDMRG calculation, the cylindrical geometry breaks the symmetry between and bonds also in the MFT, resulting in different amplitudes for different bonds. Therefore, in Figure 9 we plot the excitation spectrum of for , with larger bond anisotropy for larger . The left panels show the spectrum of and Majorana fermions, while the right panels show the minimal energies required to excite two Majorana fermions, as given by

For finite anisotropy, opens a gap at all allowed momenta, and consequently, also in . Nevertheless, the and points remain soft, as in the TM spectrum. Furthermore, additional soft modes appear at finite in the spectrum. The positions in momentum of these additional soft modes changes as a function of . In several cases, the positions of these minima are in fact similar to the positions of the soft modes we found in our analysis of the TM calculation. Compare, for example, the plots in Figure 9 with the plot of Figure 5. Finally, we note that also in the MFT, the single and two particle modes become soft at the same momenta (except at ).

5Dynamic structure factor

5.1Zero magnetic field

Next, we turn to spectral signatures of the spin liquid, which can be observed in experiments. The dynamic structure factor, which is probed in inelastic neutron scattering experiments, is defined as

We calculated using an ED method for a 24-site cluster, see Appendix C for details. Figure 10 shows for several values of . The first evident feature is the existence of a broad excitation continuum at high frequencies. In the Kitaev limit, , most of the spectral weight is found at relatively low energies. As is increased, the low frequency spectral weight at the Brillouin zone center ( point) is pushed towards higher frequencies, while a low signal remains at the , and (midway between and ) points. The difference in momentum dependence between high and low frequencies is clearly evident by integrating over different ranges of .

Figure 10: Intensity plot of the dynamic structure factor, S({\bf q},
    \omega), presented along high symmetry lines of the Brillouin zone for several values of \phi, as obtained using the ED method. Pseudo-color indicates relative intensity(blue - low, red - high). Inset: Brillouin zone with labels of symmetry points used in this figure.
Figure 10: Intensity plot of the dynamic structure factor, , presented along high symmetry lines of the Brillouin zone for several values of , as obtained using the ED method. Pseudo-color indicates relative intensity(blue - low, red - high). Inset: Brillouin zone with labels of symmetry points used in this figure.

In Figure 11 we show the relative intensity of , integrated over low and high frequency ranges. In the Kitaev limit, , is rather featureless. As is increased, , integrated over a low frequency range, shows a star shaped pattern, similar to the pattern seen in the -RuCl neutron scattering experiments at low energies. In contrast, integrating over a range of higher frequencies, shows an almost featureless momentum dependence even for finite , again, in qualitative agreement with the experiments.

Figure 11: Dynamical structure factor, S({\bf q},\omega), integrated over low (top) and high (bottom) frequencies, as obtained from exact diagonalization. The discrete set of peaks in momentum space, resulting from the finite size of the studied cluster, was broadened in order to improve the visualization of the obtained pattern. The black lines depict the boundaries of the first and second Brillouin zones.
Figure 11: Dynamical structure factor, , integrated over low (top) and high (bottom) frequencies, as obtained from exact diagonalization. The discrete set of peaks in momentum space, resulting from the finite size of the studied cluster, was broadened in order to improve the visualization of the obtained pattern. The black lines depict the boundaries of the first and second Brillouin zones.
Figure 12: Equal-time structure factor S({\bf q}), obtained using the Majorana mean-field theory.
Figure 12: Equal-time structure factor , obtained using the Majorana mean-field theory.

To calculate the dynamic spin structure factor in the context of the Majorana MFT, one must consider Z flux excitations with respect to the ground state, since each spin operator inserts a flux[22]. Technically, this requires solutions to Eqns. (Equation 5) and ( ?) which go beyond the uniform ansatz considered here. It is however still possible to approximately calculate the equal-time spin structure factor,

Since , and noting that according to the ED results, most of the dynamic structure factor signal is concentrated at low frequencies, we expect that resemble the integrated low frequency patterns of . Indeed, as seen in Figure 12, the mean-field theory reproduces the star shaped pattern seen both in experiments and in the ED calculation.

Figure 13: Dynamic structure factor for \phi/\pi=0.2, J_3=0.05, and several values of the magnetic field h_{\perp c^*}. Symmetry points labeled as in Fig. .
Figure 13: Dynamic structure factor for , , and several values of the magnetic field . Symmetry points labeled as in Fig. .
Figure 14: Magnetization and second derivative of the energy as a function of magnetic field, parallel , h_{||c^*}, and perpendicular, h_{\perp c^*}, to c^* = (1,1,1).
Figure 14: Magnetization and second derivative of the energy as a function of magnetic field, parallel , , and perpendicular, , to .

5.2 terms and finite magnetic field

Zig-zag magnetic ordering, similar to the magnetically ordered state observed in -RuCl at low temperatures, can be stabilized by adding a third neighbor Heisenberg term, to the Hamiltonian, Eq. (Equation 1). Furthermore, it is possible to suppress this ordering tendency by applying a magnetic field, . This is evident in Figure 13, which displays for , and several values of the in-plane magnetic field , perpendicular to , which corresponds to the out-of-plane direction. With the low spectral weight is increased at the M and Y points, but not at . Notice, however, that most of the zone center ( point) spectral weight is found at relatively high energies. This aspect of the spectra is similar to the case with , , shown in Figure 10. When is increased beyond , the zone center spectral weight shifts towards lower energies, and a continuum of excitations emerges.

Figure 14 shows magnetization curves as obtained with ED, for magnetic fields pointing parallel, , and perpendicular, , to . A peak at in the second derivative of the energy , indicates an apparent transition away from zig-zag order. In addition, the magnetization curves display an easy axis anisotropy, also consistent with -RuCl experiments[12]. A simple mean-field analysis qualitatively explains the easy-plane anisotropy as follows: If a ferromagnetically ordered moment is assumed as a classical Weiss field, the mean-field energy is obtained as

which is minimized, for finite , when is satisfied, i.e., when the magnetic moments are in-plane.

6Conclusions

In this work, we investigated a spin model with both the Kitaev () and symmetric-anisotropic () interactions on the honeycomb lattice using iDMRG, exact diagonalization, and Majorana mean-field theory. This model is strongly motivated by recent experiments on -RuCl, where and are likely to be the dominant exchange interactions.

We found strong numerical evidence for the existence of a quantum spin liquid for arbitrary ratio of for ferro-like Kitaev interactions in iDMRG. In particular, the entanglement entropy remains very high in this entire region while we do not see any sign of magnetic order in iDMRG computations. In contrast, we found a magnetically ordered state with very small entanglement entropy on the antiferro-like Kitaev side. Moreover, we demonstrated the existence of coherent two-dimensional multi-particle excitations using the correspondence between transfer-matrix eigenvalues and the lower boundary of multi-particle excitation spectrum. The cylinder geometry in iDMRG induces an anisotropy in bond-dependent energy, which is expected to move the locations in momentum space of low energy excitations. We show that this can indeed be seen in the transfer-matrix spectra. The existence of such two-dimensional coherence excitations without magnetic order is a very strong evidence of quantum spin liquid.

In order to make direct connection to neutron scattering experiments, we computed the dynamical structure factor for the model without and with a small third neighbor Heisenberg interaction using exact diagonalization of 24-site cluster. The on top of the interactions is shown to drive a transition to the zig-zag order in agreement with previous numerical studies. Upon introduction of starting from the ferro-like Kitaev interaction, the dynamical structure factor develops the scattering continuum with star-like intensity profile at low energies, just like what is seen in recent neutron scattering experiment. The magnetic field dependence of the dynamical structure factor is also investigated when is finite such that the ground state is the zig-zag magnetic order in zero field. There is a transition to a paramagnetic state with dominant scattering intensity at the zone center when the external magnetic field along the honeycomb plane reaches about of the largest exchange interactions, namely or . This is seen in the magnetization profile and dynamical structure factor computed in exact diagonalization. All of these features are consistent with recent experimental data.

Further we used the Majorana mean-field theory to gain analytical insight in these numerical results. For example, we computed single and two-particle excitation spectra for the model and compared these with transfer-matrix eigenvalues in iDMRG. This provides a rule of thumb for understanding the multi-particle excitation spectrum seen in transfer-matrix. The equal-time structure factor and real-space spin correlations computed in the Majorana mean-field theory are also consistent with main features in the exact diagonalization results.

Combing all these results together, we conclude that the mysterious scattering continuum seen in the neutron scattering experiment on -RuCl may come from a nearby quantum spin liquid supported by interactions. In our numerical computations, the spin liquid phases at finite show qualitatively the same behavior as the Kitaev spin liquid. We do see, however, a jump in the bond-dependent energy at some value of in iDMRG on cylinder geometry, which causes a small kink in the entanglement entropy. This can be interpreted as a “meta-nematic” transition, where the bond-anisotropy (or broken 3-fold rotation symmetry) increases abruptly. Whether such a transition would survive in 2D limit is not clear to us at present. We emphasize that the quantum spin liquid phases before and after this jump appear to be quite similar judging from the transfer-matrix spectra. Hence even if it survives in 2D limit, it would be more like a Lifshitz transition on the Fermi surface topology of the underlying quasiparticles. These questions will have to be addressed in future theoretical investigations. Further experimental data in external magnetic field would provide additional clues for the validity of the assumption that the or is a good minimal model for -RuCl.

Acknowledgements

We would like to thank Andrei Catuneanu, Yin-Chen He, Masatoshi Imada, Hae-Young Kee, Young-June Kim, Stephen Nagler, Tsuyoshi Okubo, Natalia Perkins and Roser Valenti for useful discussions. YY further thanks Mitsuaki Kawamura for his technical support. YY was supported by PRESTO, JST (JPMJPR15NF). The ED computation has partly been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo. MG and FP acknowledge support from the German Research Foundation (DFG) via SFB 1143 and Research Unit FOR 1807. GW and YBK are supported by the NSERC of Canada, Canadian Institute for Advanced Research, and the Center for Quantum Materials at the University of Toronto.


AiDMRG and transfer matrix spectrum

This appendix is devoted to the infinite Density Matrix Renormaliation (iDMRG) method and the transfer matrix spectrum. The first part exposes the geometry used in order to apply iDMRG onto a 2D lattice model. We will present additional data to characterize possible finite size effects, due to the finite circumference. The second part introduces the transfer matrix spectrum and how to make use of the ground state obtained by iDMRG.

Infinite Density Matrix Renormalization Group.

We use iDMRG [17] to study ground state properties of the K- model. Initially developed for 1D system, it has been successfully used for 2D systems by wrapping the lattice on a cylinder and mapping the cylinder to a chain. Furthermore employing translational invariance allows to study infinite cylinders[18]. Due to the cylinder geometry, one dimension of the lattice is finite and leads to a discretization of the related reciprocal vector. Thus, accessible momenta lie on lines in reciprocal space. We chose cylinder geometries, i.e. circumference and unit cell, such that those lines go through the gapless nodes of the isotropic Kitaev spin liquid, which are located at the -points of the first Brillouin zone. The results presented here and in the main text are obtained using a rhombic unit cell and a narrow cylinder with sites circumference.

We extend the discussion of the main text by considering the average energy of the and bonds for , which is presented in Figure 6. In the limit of small coupling almost no anisotropy exists, which indicates negligible finite size effects caused by the cylinder geometry. Once the -coupling increases, the anisotropy raises and reaches near . Using wider cylinders with reduces the anisotropy only slightly. This suggests, that a -like interaction is highly sensitive to finite size and a cylinder geometry and possibly causing the gapless Dirac nodes to shift away from the K-points in reciprocal space.

Transfer Matrix Spectrum.

The transfer matrix (TM) of an infinite matrix product state (iMPS) contains full information about the static correlations along the wavefunction encoded as iMPS[20]. Intuitively, the TM translates the iMPS by one lattice vector along the chain in 1D or cylinder in 2D. The static correlations are related to the spectral gap[32], e.g. for . This statement has been extended by Zauner et al. [20] to also include momentum, such that the lengthscale of the decay of static correlations with a certain momentum gives an upper bound on the spectral gap with close to . Thus the transfer matrix spectrum, a ground state property, provides certain information about the dispersion, such as where the minimal gap is situated or the position of a potential Dirac cone.

Let be the eigenvalues of the transfer matrix with the ordering . By definition the dominant eigenvalue is . Generally, are complex and can be decomposed as . The angle is connected to the momentum along the chain or cylinder. Exploiting the rotational invariance of the Hamiltonian on the cylinder geometry yields the transverse momentum as will be explained now. In the following, we require the iMPS to be in canonical form[33]. A translation with a lattice vector along the circumference keeps the Hamiltonian invariant and as such can be treated as a regular quantum number. We extract by computing the dominant eigenvector of the mixed transfer matrix constructed out of the ground state iMPS and a iMPS with the translation applied, see also Figure 17b). We like to remark, that the translation along the circumference is simply given by a permutation of sites within a ring. If the iMPS is sufficiently converged and the applied translation is indeed a symmetry, then the dominant eigenvalue of the mixed TM is . Its eigenvector has a diagonal form with eigenvalues and being discretized in steps , where is the number of unit cells around the cylinder. If some Schmidt values are degenerate, the diagonal form becomes block diagonal with blocks for each set of degenerate Schmidt values. Each block can be diagonalized separately by a unitary transformation which is then applied to the non-translated iMPS. The momentum quantum number are associated with the entries along an bond leg in the same way as Schmidt values are. The TM connects states and with corresponding and , hence the transverse momentum is given by . The label of can be read off from its eigenvector due to the fact, that has only non-zero entries with the same change of the quantum number .

Figure 15: Schematic representation of a) the regular and b) the mixed transfer matrix \mathcal{T} with translation T applied. c) Dominant eigenvector \tilde \Lambda_0 of \mathcal{T}^T determines the q quantum numbers associated with each bond leg.
Figure 15: Schematic representation of a) the regular and b) the mixed transfer matrix with translation applied. c) Dominant eigenvector of determines the quantum numbers associated with each bond leg.


Figure 16: Schematic representation of a) the regular and b) the mixed transfer matrix \mathcal{T} with translation T applied. c) Dominant eigenvector \tilde \Lambda_0 of \mathcal{T}^T determines the q quantum numbers associated with each bond leg.
Figure 16: Schematic representation of a) the regular and b) the mixed transfer matrix with translation applied. c) Dominant eigenvector of determines the quantum numbers associated with each bond leg.


Figure 17: Schematic representation of a) the regular and b) the mixed transfer matrix \mathcal{T} with translation T applied. c) Dominant eigenvector \tilde \Lambda_0 of \mathcal{T}^T determines the q quantum numbers associated with each bond leg.
Figure 17: Schematic representation of a) the regular and b) the mixed transfer matrix with translation applied. c) Dominant eigenvector of determines the quantum numbers associated with each bond leg.

BMajorana MFT details

Even before solving equations (Equation 5) and ( ?), it is important to characterize the behavior of the Majorana fermions under any configuration of and . The operators describe fermions which move about the whole lattice with hopping amplitudes given by . Similarly, the operators describe fermion hopping with amplitudes . Most importantly, the structure of , given in Eq. (Equation 3), separates the fermions into three independent, uncorrelated sectors. Each sector is associated with one of the three sublattices of hexagons. Within each sector, the hopping amplitude is around a hexagon of the corresponding sublattice, and between neighboring hexagons. When , fermions are bound to a bond, while they are bound to a hexagon when . This behavior echoes the macroscopic degeneracies of the parent classical models: in the classical Kitaev model there are a macroscopic number of degenerate ground states which are related to each other by reversing the sign of a single spin component for two spins on the same bond. The classical -model has a similar macroscopic degeneracy, obtained by reversing the signs of one component of each of the six spins on one hexagon[34]. Furthermore, returning to the uniform mean-field solution, the -fermion band structure is similar to the band structure obtained in a Luttinger-Tisza study of the corresponding classical models, since both are determined by the exchange matrix . In the classical case, the flat bands for or indicate the existence of a macroscopic number of ground states, mentioned above.

Figure 18: Spin-spin correlations in real space. The open black circles indicate site i=0 and the filled circles indicate i. The circle size and color represent the magnitude and value of the corresponding correlation \braket{S_i^xS_0^x}. (a) Mean-field results for \phi/\pi = 0.4. (b) ED results for \phi/\pi = 0.5. In both cases one finds same spin component correlations only among a subset of all sites.
Figure 18: Spin-spin correlations in real space. The open black circles indicate site and the filled circles indicate . The circle size and color represent the magnitude and value of the corresponding correlation . (a) Mean-field results for . (b) ED results for . In both cases one finds same spin component correlations only among a subset of all sites.

Spin correlations in real space, shown in Figure 18a, reveal a pattern which manifests the separation of fermions into independent sectors: on site is correlated only with a certain subset of ’s on other sites . Similar patterns are also obtained with ED, see Figure 18b although some differences should be pointed out. According to the MFT, at (), there are non-zero spin-spin correlation only within the same hexagons since the -fermions are localized. Thus, both finite and are requried to get longer range correlations. Furthermore, the correlation in Figure 18b don’t decat very fast since they were obtained using ED on a small cluster with periodic boundary condition, whereas the correlations in Figure 18a where calculating assuming an infinite system. Finally, even among the the subset of sites which have same spin component correlations, MFT gives nonzero static correlations only between opposite sublattice sites.

CDynamical structure factors and the Krylov subspace method

The dynamical structure factor may be written as

where

and is the Fourier transform of the spin operators defined as,

The dynamical structure factors are conventionally calculated by using the Lanczos algorithm initialized with an excited state [35]. However, naive implementations of the Lanczos algorithm suffer uncontrolled truncation errors. In contrast, modern Krylov subspace methods, which extract essence from the Lanczos algorithm, offer controlled convergence.

Here, we solve a linear equation by employing a conjugate gradient (CG) method, instead of explicitly calculating the resolvent of in Eq.(Equation 9). The CG methods find the solution in a Krylov subspace, as follows. First, by introducing the following two vectors,

we rewrite as

To obtain the unknown vector , we solve the following linear equation,

When the linear dimension of the matrix , , is too large to store the whole matrix in the memory, the linear equation is solved iteratively, for example, by using the CG methods. At th iteration, the conjugate gradient algorithm initialized with finds an approximate solution within a -dimensional Krylov subspace . At each steps, the CG-type algorithms search the approximate solution to minimize the 2-norm of the residual vector,

We note that one needs to solve Eq.(Equation 10) essentially once at a fixed complex number to obtain whole spectrum . Due to the shift invariance of the Krylov subspace [36], namely, for any complex number , we can obtain from with a numerical complexity of [36]. The Krylov subspace methods utilizing the shift invariance are called the shifted Krylov subspace methods.

For the calculations of , we employ the shifted biconjugate gradient (BiCG) method implemented in a numerical library for the shifted Krylov subspace method [37]. The condition for truncating the shifted BiCG iteration is set for the following calculations with . The number of the iteration steps for satisfying the condition is typically of the order of ten thousand.

References

  1. bibitemNoStop
  2. bibitemNoStop
  3. bibitemNoStop
  4. in (, ) pp. bibitemNoStop
  5. bibitemNoStop
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
21452
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description