Characterizing spin-one Kitaev quantum spin liquids

Material realizations of the bond-dependent Kitaev interactions with $S$=1/2 local moments have vitalized the research in quantum spin liquids. Recently, it has been proposed that higher-spin analogues of the Kitaev interactions may also occur in a number of materials with strong spin-orbit coupling. In contrast to the celebrated $S$=1/2 Kitaev model on the honeycomb lattice, the higher-spin Kitaev models are not exactly solvable. Hence, the existence of quantum spin liquids in these systems remains an outstanding question. In this work, we use the density matrix renormalization group (DMRG) methods to numerically investigate the $S$=1 Kitaev model with both ferromagnetic (FM) and antiferromagnetic (AFM) interactions. Using results on a cylindrical geometry with various circumferences, we conclude that the ground state of the $S$=1 Kitaev model is a quantum spin liquid with a $\mathbb{Z}_2$ gauge structure. We also put a bound on the excitation gap, which turns out to be quite small. The magnetic field responses for the FM and AFM models are similar to those of the $S$=1/2 counterparts. In particular, in the AFM $S$=1 model, a gapless quantum liquid state emerges in an intermediate window of magnetic field strength, before the system enters a trivial polarized state.


I. INTRODUCTION
The discovery of several candidate materials for the S=1/2 Kitaev quantum spin liquids [1] has been a driving force in topological quantum materials research. In this regard, the exactly solvable quantum spin liquid state of the S=1/2 Kitaev model [2] has not only settled the long theoretical debates on the existence of quantum spin liquids, but also offered practical platforms for materials realizations. Recently, a number of candidate materials were proposed for the higher-spin S=1 Kitaev systems [3], anticipating possible quantum spin liquids. The higher-spin Kitaev models are not exactly solvable, and hence the nature of quantum ground states is currently an important subject of theoretical investigation.
Some basic properties of the S=1 Kitaev model on the honeycomb lattice are known [4]. In analogy to the S=1/2 model, one can define a plaquette flux operator W p on each hexagon, which commutes with the Hamiltonian. These plaquette operators are written in terms of the π-spin-rotation unitary operators for an arbitrary spin S quantum number [4]. Hence, each eigenstate can be labelled by the eigenvalues of these flux operators. This constant of motion can be used to show that spinspin correlations exist only for nearest neighbours, however this by itself is not enough to allow for exact solutions. Exact diagonalization (ED) studies up to 24-site clusters [5] concluded that the ground state might be a gapless quantum spin liquid. On the other hand, a recent tensor network construction of the variational wavefunction [6] suggests that the ground state is a gapped Z 2 spin liquid with abelian quasiparticles. Since each numerical approach has its own limitations, it is important to synthesize the efforts from different numerical and analytical approaches for the ultimate understanding.
In this work, we study the S=1 Kitaev model on the honeycomb lattice using DMRG [7][8][9] on a cylindrical FIG. 1. Different geometries for the honeycomb Kitaev model. Blue, green, and red lines represent the x-bonds, y-bonds and z-bonds respectively. Panel A is the two-leg ladder (Ly = 2), which is equivalent to a square ladder. Panel B shows the three-leg ladder (Ly = 3) with periodic boundary conditions along one of the axes. Red dashed lines represent the periodic links, and roman letters represent the periodic z-bonds. geometry with various circumference lengths. We start with two-leg ladder systems (or L y = 2), where we study system sizes up to 250 sites (L x = 125). The ground state has a uniform flux, for which W p = 1 for any hexagonal plaquette. We demonstrate that the spin-flip operator at a given site generates two adjacent "vortex" plaque-arXiv:2001.06000v1 [cond-mat.str-el] 16 Jan 2020 ttes with W p = −1, just like in the S=1/2 case. From the two-fold degeneracy of entanglement spectrum (ES), we conclude that the ground state of the two-leg ladder system is a symmetry-protected topological (SPT) phase, with a two-fold degenerate ground state. This result is similar to the two-leg ladder system of the S=1/2 model [15], albeit the degeneracy structure of the ES in the S=1 model is different from that of the S=1/2 model. Given that the S = 1 model naturally offers an AFM Kitaev exchange interaction unlike the J eff = 1/2 FM Kitaev model, it is worthwhile to investigate the phase diagram under the magnetic field. We apply fields perpendicular to the honeycomb plane (parallel to the [111] direction). In the S=1/2 model [10], it is known that the magnetic field phase diagram of the two-leg ladder system is very similar to that of clusters with wider circumferences [11]. Examining the magnetic field responses of the FM and AFM Kitaev couplings, we find that the magnetic field dependence is surprisingly similar to the S=1/2 case [12]. For example, an intermediate phase exists for the AFM model in a window of magnetic field strengths, right before the system enters a polarized state, while for the FM model there is a direct transition to the polarized state at a much lower critical field. On the other hand, in comparison to the S=1/2 model phase diagram, the zero-field ground state of the AFM model is much more robust against the magnetic field than that of the FM case.
We then consider a three-leg (or L y = 3) system with periodic boundary conditions along the circumference, with cluster sizes of up to 144 sites (L x = 48). The ES no longer shows any degeneracy. Furthermore, the ground state energy of the FM and AFM models is exactly the same, just like for the S=1/2 model. We define the Wilson loop operator along the circumference direction W , which commutes with the Hamiltonian. We further investigate the ground state of the AFM model, and its two lowest energy excited states. The first three lowest energy states are in the W = +1 sector. The lowest excitation energy density (the difference, per site, between the first excited state energy and the ground state) for the largest system at hand, L x = 24, is ∆e = 7×10 −4 K, where K is strength of the Kitaev interaction. This is our upper bound on the excitation gap of the S=1 AFM Kitaev model, which may be considered to be quite small. We also study how the ground state topological properties change with increasing the number of legs. We use systems of up to 48 sites with L y ≤ 6, and find that the ground states of L y = 4 and 6 clusters are in the W = −1 sector, while for L y = 3 and 5 clusters it is in the W = +1 sector . Therefore it appears that there is an even-odd effect, which suggests that the ground states in W = +1 and W = −1 sectors may become degenerate in the thermodynamic limit. This would be consistent with the two degenerate ground states on the cylinder for Z 2 spin liquid states.
The effect of external magnetic fields is studied for the three-leg cylinder consisting of 24-sites. The overall be-haviour is similar to the results of the two-leg ladder system. That is, the critical field for the AFM model is much larger than that of the FM case. Due to slow convergence, however, we can only see the first transition to the intermediate state in the AFM model for the L y = 3 system. The transition is also well apparent in other observables such as the entanglement entropy (EE). We investigate the excitation energy gap as a function of an external magnetic field and find that it vanishes as one approaches the aforementioned transition. This is consistent with the picture of a gapless intermediate phase, just like for the S=1/2 AFM model subject to magnetic field [13,14]. The general overall picture emerging from these studies is that the ground state of the S=1 Kitaev model is a quantum spin liquid with a Z 2 gauge structure, and that the response to external magnetic fields is very similar to the S=1/2 case. Our current numerical data may be consistent with a quantum spin liquid with a small excitation gap. This, however, does not preclude a gapless spin liquid in the thermodynamic limit.
The rest of the paper is organized as follows. In Section II, we introduce the model and explain basic symmetry properties. Here we also discuss the cylindrical geometry we use throughout this manuscript. In Section III, we briefly explain the DMRG details. We discuss the results of the two-leg ladder systems or L y = 2 in Section IV. In Section V, we present the results of the three-leg system L y = 3. Here we also briefly discuss more results on systems with up to six-leg ladders (L y = 6). In Section VI, we discuss the implications of our results.

II. MODEL
The Kitaev model Hamiltonian is given by where S γ i is the γ component of an S = 1 spin at site i on a honeycomb lattice, i, j γ is two nearest-neighbour sites along an γ bond (γ = x, y, z). h is a uniform magnetic field, where we discuss a field along the [111] direction.

A. Symmetries
The pure (h = 0) Kitaev model on a honeycomb lattice, which is defined in Eq. 1, has a constant of motion defined on a plaquette where P j is the j-th plaquette consisting of six sites on a single hexagon, and α represents the protruding bond along P. circumference of the cylinder, which commutes with the Hamiltonian where y-loop is a closed loop around the circumference in the periodic direction. For the S = 1/2 Kitaev model, it can be shown that spin operators S α i acting on any eigenstate of the Hamiltonian, would lead to a π-flux insertion, or a sign-flip of the associated plaquette operator W j p , where j corresponds to two adjacent plaquettes which include site i, and share an α-bond. Since W j p commutes with the Hamiltonian, it makes the spin-spin correlation function short-ranged, such that it is non-zero for nearest neighbours and exactly vanishes for further neighbours, This property remains the same for the S = 1 Kitaev model (as was previously semiclassically proven for higher spin [4]).

B. Geometries
In this study we mainly focus on two geometries, one of a two-leg ladder and one of a three-leg ladder (see Fig. 1). It was previously shown [11] that the two-leg ladder's restricted geometry is able to capture the phase transitions for the S=1/2 model fairly well. Later in [15], it was demonstrated that this simplified geometry's phase diagram in the Kitaev-Heisenberg plane is very similar to the 2D phase diagram on the honeycomb lattice. Hence, we expect that any transition to be found for the S = 1 ladder would also appear in the 2D phase diagram albeit with different critical parameters (field strength, etc.).
An additional advantage is that it is much easier, numerically, to access very large systems due to shorter-range interactions.
The three-leg ladder is the minimal geometry which allows the Wilson loop operator along the circumference W . In addition, as discussed in Ref. [16], it allows probing the high-symmetry K-points in the Brillouin zone, which host Dirac fermions in the spin-half case.

III. RESULTS
We study the spin-one Kitaev model on a honeycomb lattice, using the DMRG [7][8][9] technique. We use cylinders with open boundaries conditions, with up to 250 sites for the two-leg geometry, or L x = 125, L y = 2, (see Fig. 1 B), up to 144 sites for the three-leg geometry, or L x = 48, L y = 3, (see Fig. 1 A). We retain 7200 states in the reduced density matrix, with no symmetries kept, and we found that 45 sweeps were sufficient for good convergence. The DMRG relative truncation error was less than 10 −9 .

A. General ground state properties
We consider both the FM and AFM Kitaev couplings (K = ∓1 in Eq. 1). Note that the AFM couplings were shown to host a wider Kitaev spin liquid phase in the parameter space of the Kitaev-Heisenberg model [3]. Moreover, one should note that under addition of a sufficiently small nearest-neighbour Heisenberg interaction, the following results are still valid, and that the phase diagrams presented here remain qualitatively the same.
Wilson loop operators: W -We examine ladders of up to L y = 6. One can define the W operator (see Eq. 3) on a geometry consisting of three legs and above. We find an even-odd effect for different topological sectors. For the odd-number-leg systems we find the ground state to be in the W = 1 sector, while for the even-number-leg clusters the ground state is in the W = −1 sector.
Short range correlation effects -The spin S = 1 model hosts very short spin-spin correlations. For the pure Kitaev limit, h = 0 point, spin-spin correlations are non-zero only along the specific bonds, i.e. S α i S β j ∝ δ αβ δ i,j α . This is a direct consequence of the role of spin operators acting on eigenstates of the Hamiltonianthey thread π-flux into two adjacent plaquettes, i.e. they change the local plaquette expectation value by a factor of −1. However unlike the S = 1/2 case, where the flux insertion operators are Pauli matrices, the spin S = 1 operators alter the normalization of the state. Hence the resulting state after the spin flip should be written as the following state |ψ , where two of the adjacent plaquettes to site i, sharing an x-bond, gain the additional π-flux: eigenstate). Moreover, at finite field strengths, h, although W p is not a conserved quantity anymore, applying the spin operators still flips the sign of the two adjacent plaquttes. Furthermore, W p does not commute with the Hamiltonian, and longer-ranged correlations appear and the property that only bond-dependent correlations are present and non-zero is lost, i.e. S α i S α j = 0 for any two sites (i, j).

A. Ground state properties
The two-leg ladder, with boundary conditions depicted in Fig. 1 (B), shows doubly degenerate ground states in the thermodynamic limit (the energy difference between the ground state and the first excited state with L x = 125 is ≈ 2 · 10 −14 ). For both FM and AFM coupling the ground state energies are the same, and the two degenerate states have uniform flux, W p = 1. Interestingly, for the AFM coupling at fields below h c1 (the lower critical field marking the intermediate phase as discussed later), the degeneracy is still present, confirming the topological nature of the ground state (due to the symmetry protected topological, SPT, phase which is discussed below). At fields h < 0.2 < h c1 , the energy difference between the two lowest energy states is ≈ 10 −12 . At h = 0.2 the gap density is ≈ 3·10 −4 , which can still be attributed to finite size effects.
Entanglement spectrum and entanglement entropy -An SPT phase can be characterized by double degeneracy of its ES [17]. It is no surprise that for the two-leg ladder geometry [18] we find that the ES is doubly degenerate. In fact, at h = 0 the degeneracy pattern of the density matrix eigenvalues, λ i is 2 − 4 − 2 for the entire ES (see Fig. S1). This aforementioned degeneracy breaks in the presence of magnetic field. In addition, the two transitions to the intermediate phase, and to the polarized phase are revealed by examining the EE as seen in Fig. 5. At the pure Kitaev limit, h = 0, the entropy is high, it then jumps to a lower value in the presence of a weak magnetic field, indicating a phase transition from the Kitaev limit. As one further increases the field, the EE rises until it jumps once more at h = h c1 , indicating the intermediate phase, then there exists another jump at h = h c2 . After increasing the field beyond h = h c2 , the EE begins to drop, an indication of an order that is being built, which is the polarized phase.  Fig. 1 A). The magnetic field, whose magnitude is h, is uniform, and parallel to the [111] direction. Panel A shows the ground state energy density, panel B is the magnetization density, panel C is the uniform magnetic susceptibility. Panel D is the total spin magnitude, and (E) is the plaquette operator's, Wp, expectation value. Panel F shows the difference in energy between the first two excited states and the ground state. This quantity suggests that, similarly to the spin S = 1/2 case, the specturm collapses at the phase transition to the intermediate phase. This transition is also captured by the uniform magnetic susceptibility. Wp does not show any distinct features near the transition.

B. Magnetic phase diagram
The Magnetic phase diagram obtained for the two-leg ladder, is summarized in Fig. 3. For both FM and AFM couplings, we plot the ground state energy, the magnetization density, and the uniform magnetic susceptibility, as a function of magnetic field parallel to the [111] direction (panels A-C respectively of Fig. 3). For the FM there appears to be a single transition, at weak fields, to a polarized phase. However, for the AFM coupling, one notices two kinks around h c1 ≈ 0.34 and h c2 ≈ 0.48 in the magnetization, and two peaks at the same positions for the magnetic susceptibility. While the magnetization density parallel to [111] is increasing as expected, the perpendicular (in plane) magnetization remains negligible regardless of the strength of the magnetic field applied. The spin size, |S total | (depicted in Fig. 3 panel D) is extracted from the ground state's expectation value S 2 total . It displays a two-kink structure at the same critical fields. Unlike the two previous observables, from the plaquette operator, W p , whose expectation value is given in Fig. 3 E, one cannot identify the exact location of the phase transitions. Moreover, since in the presence of magnetic field it is no longer a constant of motion, the decrease in its value goes hand in hand with the increase of longer-ranged spin-spin correlations.
The two peaks allow us to divide the phase diagram into three regions. For h < h c1 we get a distinct phase where the magnetization starts to build up, and where the susceptibility and total spin are increasing non-monotonically. By solving the entire spectrum for small clusters using ED (see SM for more details), we see that the distance between eigenenergies begin to shrink, i.e. the density of states at low energies is increasing, and, finite size-gaps begin to decrease. Once the spectrum collapses at h c1 , a new phase appears for fields h c1 < h < h c2 . This phase is characterized by the peaks appearing in the magnetic susceptibility. This intermediate phase, which shows the least convergence, and requires much larger bond-dimensions (> 7200 kept states), could possibly have a diverging correlation length, and hence hinting that it might be a gapless phase -consistent with the spectrum collapse. Finally, a polarized phase appears for fields h > h c2 , where a linear increase in the magnetization and the total spin is seen.   (Fig. S1) for the entanglement spectrum.

A. Ground state properties
We start by examining the three lowest eigenstates of the three-leg ladder geometry. The energies of these states are portrayed in Fig. 2. One can see that as system size increases, the gaps between these states decrease.
We find that the ground state and the next two eigenstates are always in the W = 1 sector. As shown in Fig. 2, the energies of these three states become very close as the system size increases. Unfortunately, based on the finite size date obtained here, we cannot definitely determine whether the system has a degenerate ground state followed by an excitation gap or a gapless spectrum in the thermodynamic limit.
Similar to previous ED results we find that as one increases the number of legs the ground state energy density actually increases, e.g., for the two-leg ladder we find e 0 ≡ E0 N = −0.672 K, for the three-leg ladder e 0 = −0.644 K.

B. Magnetic phase diagram
The magnetic field dependence of the ground state of the three-leg system is depicted in Fig. 4. We plot the ground state energy density, the magnetization density, the uniform magnetic susceptibility, the total spin magni-tude, and the plaquette operator expectation value. We perform this analysis for both the FM and AFM interaction couplings. For the ferromagnetic case, the field undermines the Kitaev interaction and easily polarizes the system (as can be seen in panels B and D), therefore we conclude that the transition to the polarized phase occurs at weak fields. However, for the AFM coupling, the field has to compete with the staggering nature induced by the Kitaev interaction. Therefore, the magnetization builds up much slower, and the total spin is increasing slowly as a function of the applied field (panel D). For the AFM case, the transition to the intermediate phase is seen in panel C by the wiggling of the uniform magnetic susceptibility and in panel F where we show that first two excited states energies decrease towards the ground state energy, which might indicate a spectrum collapse (similar to the two-leg ladder, see SM for ED spectrum).
Entanglement spectrum and entanglement entropy -The EE for the three-leg ladder as a function is summarized in Fig. 5. Unlike for the two-leg ladder system, there is no jump in the ES at the Kitaev limit (h = 0), implying there are no symmetry protected features in the EE unlike the two-leg geometry. The latter claim is also supported by the ES in Fig. S2, which shows no degeneracy. In addition, as seen in Fig. 5, the EE is monotonically building up towards the transition to the intermediate phase at h = 0.22.
We expect that this qualitative picture of the two transitions will not change in systems with larger circumference. However, based on studies of smaller clusters, the critical field values are expected to be somewhat affected by finite size effects.

VI. DISCUSSION
In the current work we used DMRG to investigate the nature of the ground states for the S=1 Kitaev model with both FM and AFM exchange interactions. We presented these results via systematically studying the evolution of the ground states as a function of the circumference size L y of the cylinder, and applied magnetic field. We found a number of numerical evidence suggesting that the ground states of these models for L y ≥ 3 are quantum spin liquids.
First, we have identified an SPT for the S = 1 twoleg ladder geometry via the two-fold degeneracy of the ground state and the degeneracy of the ES. As one introduces magnetic field for the AFM model, three phases are found: a disordered, highly-entangled phase at weak fields, a gapless intermediate phase and a polarized phase. For the FM couplings, a direct transition to the polarized phase is found at a weaker field. We find a similar magnetic phase diagram for the three-leg ladder, and argue that due to the similarity with the two-leg ladder system, both would share the same magnetic responses.
For the three-leg ladder, we find an upper bound on the excitation energy of the AFM model -∆e = 7 × 10 −4 K, which may suggest a quantum spin liquid with a small excitation gap or a gapless spin liquid in the thermodynamic limit. However, given our data, we cannot distinguish between a gapless spin liquid or an existence of a small gap.
We find great similarity between the S = 1/2 and S = 1 models: (i) The ground-state is two-fold degenerate for the two-leg ladder. (ii) Both models share a similar response to applied magnetic field. In the AFM case, we obtained a phase diagram which is separated into three distinct regions (although the critical fields values differ by a factor of ∼ 1.5 between S = 1/2 and S = 1 models). The critical field strengths marking the phase transition are accompanied by a spectrum collapse, which suggests a gapless disordered state. (iii) The ground state energies of the FM and AFM couplings are the same. (iv) Both models exhibit extremely short-ranged correlations. The correlations become longer-ranged with applied field, together with a gradual drop in the magnitude of the (local) plaquette operators.
Still, there are differences between the two models. The even-odd effect we find for the S = 1 model with respect to the ground state's W sectors does not occur for the spin S = 1/2 model. Another difference lies in the ES structure, while the S = 1/2 two-leg ladder has fourfold degeneracy in its ES, the S = 1 ES has a 2 − 4 − 2 structure [19]. This could be due to a different symmetry protecting the SPT. Furthermore, note that for the soluble S = 1/2 case the eigenvalues of the Wilson operators correspond to Z 2 fluxes, and the ground state is in the W = ±1 sector depending on the boundary conditions (periodic versus anti-periodic boundary conditions for the Majorana fermions [16]).The same flux sectors are present for the S = 1 model, as was shown in this paper, however, it is still unclear what the nature of the elementary excitations in this model is.
Using the intuition previously obtained from spin S = 1/2 calculations, we would like to point out that the twoleg ladder results may be used as a guide at a qualitative level to the two-dimensional limit. This is also seen by comparing the three-leg ladder phase diagram with that of the two-leg ladder in the presence of a magnetic field. We expect that for the AFM S = 1 Kitaev model, in the thermodynamic limit, there would be three distinct phases as well, with a gapless intermediate phase, as suggested by our results. For the multi-leg ladders, the evenodd effect of the ground state's W subspace suggests that in the thermodynamic limit, both W = ±1 sectors may be degenerate, hence it would restore a two-fold degeneracy of the ground-state on a cylinder.
Based upon all the physical quantities we discussed above it is highly suggestive that the ground state of S = 1 Kitaev model is a quantum spin liquid with a Z 2 gauge structure.
An open question remains whether a gap opens in the presence of a magnetic field. As was pointed out, at zero field, the gap could be vanishing or small. However, by viewing the spectrum in the presence of a field (Fig. 4 F), one could conclude that the finite size gap is decreasing even further. If this is correct, it could mean that in the thermodynamic limit, the entire Kitaev phase is gapless. However, a different option is that similarly to the S = 1/2 model, a small gap opens as the magnetic field is introduced, while later vanishes at finite magnetic field strength. Our current data cannot distinguish between these two scenarios, and this question is left for future study.
Note added: during the completion of this manuscript, we became aware of a similar analysis [20] on the spin S = 1 Kitaev model. We thank Donna N. Sheng for the correspondence and discussion about theirs and ours results. model [22][23][24][25]. These ED and TPQ studies show the characteristic two-peak structure of the specific heat for strong J, the excitation energy, and the overall phase diagram as a function of the inter-layer coupling J. Here, we would like to target the S = 1 limit of this model. We argue that the phase transition from two decoupled (or weakly coupled) layers occurs at finite values of J. As can be seen in Fig. S3, we calculate the EE and ES of the coupled-S = 1/2 bilayer two-leg ladder as a function of the ferromagnetic inter-layer coupling, J. One can see that for finite values of J/K 10, we converge  to the same state as we would get for the pure S = 1 Hamiltonian (comparing the values reported in Fig. 5 and Fig. S1). This suggests that the phase transition from two decoupled, or even weakly coupled S = 1/2 chains occurs at finite coupling strength.
FIG. S4. ED spectrum for 18-site two-leg ladder, as a function of uniform magnetic field h. A spectrum collapse can be seen in the vicinity of h = 0.3 K.
FIG. S5. ED spectrum for various two-leg ladder with sizes N and periodic boundary conditions (PBC). The degeneracy of the lowest eigenstates is marked by the numbers next to the markers. Overall a trend towards two-fold degeneracy of the ground state is seen starting from N ≥ 12.

Eq. 3 is invalid for this geometry)
and similarly where XZ-path runs along x and z-bonds (and similarly for YZ-path), as depicted in Fig. S6. For the two-leg ladder, the degenerate ground states are distinguishable by O x , O y , and W z . In a finite system (where finite gap still exists between these two states), all the operators mentioned are +1 for the lowest eigenstate, and −1 for the next eigenstate.
It is interesting to note another commuting operator, which is the open string operator along x-direction (the axis of the cylinder) where x-string is an open string along the x-direction, which includes every spin from one end of the finite cylinder to the other end. In a torus geometry (periodic boundary conditions along x-direction), W z will become the Wilson loop operator for a closed loop along the xdirection. However, its significance for the cylinder could be further studied. We find that for the three-leg geometry the ground state is always in the W = 1 sector and W z = 1. For all the cluster sizes we examined, the first excited state also has W = 1 and W z = 1. The next excited state, however, shows W = 1, W z = −1.

YZ path
FIG. S6. Two-leg ladder (Ly = 2) geometry, which is equivalent to a square ladder. "YZ path" is used to define the Ox string operator in Eq. E1.