Magnetic Field Induced Spin Liquids in S=1 Kitaev Honeycomb Model

We investigate the ground state properties of the spin S=1 Kitaev honeycomb model under a magnetic field based on the density matrix renormalization group calculation. With the time reversal symmetry breaking due to the magnetic field, a gapped Kitaev spin liquid is identified for both ferromagnetic (FM) and antiferromagnetic (AFM) Kitaev couplings. The topological nature of such Kitaev spin liquid is manifested by the nearly quantized Wilson loop, degeneracy in the entanglement spectra and existence of edge modes. While the FM Kitaev spin liquid is destroyed by a weaker magnetic field $H_*^\text{FM}$, the AFM one demonstrates a robustness up to an order of magnitude larger critical field $H_*^\text{AFM}$. Moreover, an intermediate nonmagnetic phase appears only for the AFM case at larger fields, $H_*^\text{AFM}<H<H_{**}^\text{AFM}$, before the transition to a high-field polarized paramagnet. The stability of the Kitaev spin liquid against the Heisenberg interactions is also examined. Our findings may further inspire the investigation of recently proposed S=1 Kitaev materials.

Remarkably, the recent observation of half-integer quantized thermal Hall conductivity in α-RuCl 3 under magnetic field offers a smoking-gun evidence of the fractionalized topological state [29], leading to substantial activities in examining the magnetic field induced QSLs in the Kitaev models and materials . Although numerous efforts have been devoted to S=1/2 case, the magnetic field induced phases in the S=1 Kitaev model and the corresponding experiments are of equal interest and importance. Different from S=1/2 Kitaev model, the nature of the ground state for the pure S=1 Kitaev model is still controversial [52][53][54][55][56][57]. Recently, a mechanism for realizing the S=1 Kitaev interactions by considering both SOC and strong Hund's couplings in layered transition metal oxides A 3 Ni 2 XO 6 (A=Li, Na and X=Bi, Sb) has been proposed [58], providing a platform to explore the S=1 Here, γ = x, y, z represent the three nearest-neighbor links i, j of the honeycomb lattice, S γ denotes effective spin S=1 degrees of freedom sitting on each site and interacting via Kitaev coupling K γ . The second term is the Zeeman term with a uniform magnetic field, and the third term represents the Heisenberg interactions with coupling J. We consider isotropic coupling K γ ≡ K and set |K| as the unit of energy, K = +1 (K = −1) corresponds to AFM (FM) Kitaev model. Below we will identify the magnetic field induced phases at J = 0 and examine the stability of KSL against Heisenberg perturbations later. We use DMRG to study the Hamiltonian (1) on a honeycomb lattice spanned by length vectors L x = L x e x and L y = L y e y , where e x = (1, 0) and e y = (1/2, √ 3/2) are two primitive vectors. As each unit cell contains two sites, and the total number of sites N = L x × L y × 2. In the present calculation, we choose both the cylinder and torus geometries and keep up to 1600∼2800 states for good convergence (depending on the nature of states and system sizes). We keep the truncation error to be smaller than 10 −6 in most cases.
We shall identify a gapped Kitaev spin liquid (KSL) phase for both antiferromagnetic (AFM) and ferromagnetic (FM) Kiatev couplings in weak magnetic fields. Interestingly while the FM gapped KSL is fragile against the magnetic field, the AFM gapped KSL is extremely robust up to an order of magnitude larger critical field. We also find the topological nature of this gapped KSL, including the Z 2 gauge structure, the edge modes and a robust multiple degeneracy in the entanglement spectra on ladder systems. Moreover, we discover an intermediate phase (IP) only in the AFM model before a second transition to the high-field polarized paramagnet. Our probe indicates that the IP may be a nonmagnetic critical regime with gapless excitations. The phase diagram is depicted in Fig. 1 (a). Finally we show that these magnetic field induced phases eventually give way to magnetic orders under sufficient strong Heisenberg-type perturbations.
Phase Diagram.-We begin with establishing the phase diagram of the model Hamiltonian (1) as a function of magnetic field strength [see Fig. 1 (a)]. We obtain the ground state for a given magnetic field, and the phase boundaries are determined by consistent evidence from the measurement results of magnetization M , magnetic susceptibility χ and the Wilson loop W y . Figures 1 (b-c) show the magnetization M and magnetic susceptibility χ = ∂M/∂H for systems with both toroidal and cylindrical boundary conditions. At the FM Kitaev coupling side [see the left-hand side in Fig. 1], we find a strong response of the system to magnetic field, the single kink in magnetization M [see Fig. 1 (b)] and single peak of χ [see Fig. 1 (c)] demonstrate a single phase transition to the highfield partially-polarized paramagnet at H FM * ∼ 0.01, which smoothly connects to the fully polarized phase in infinite field limit. In contrast, in the case with AFM Kitaev coupling [see the right-hand side in Fig. 1], the system exhibits a weak response to an external field with a critical field H AFM * ∼ 0.24, which is larger than H FM * by more than an order of magnitude [see Fig. 1 (c)]. Moreover, the two-peak structure in the magnetic susceptibility χ suggests that there is an intermediate phase at H AFM * < H < H AFM * * before a second transition to partial polarized paramagnet beyond H AFM * * . We also note that the ground state of S=1 Kitaev model at zero field can be characterized by a conserved Z 2 quantum number determined by the Wilson loop operatorsŴ l , which is generalized from S=1/2 model and can be defined i ] in every closed loop on the lattice, where α i represents the normal direction of site i [52]. For example,Ŵ y denotes the Wilson loop operator with the loop winding once around the cylinder which only cover γ = x, z links [see the inset in Fig. 1 (d)]. It is straightforward to verify that these Wilson loop operators commute with each other and also commute with the Hamiltonian of the pure Kitaev model. SinceŴ l squares to identityŴ 2 l = 1, the eigenvalues ofŴ l are W l = ±1 at H = 0, which is associated with flux excitations like in the spin-1/2 Kitaev model. At H = 0, the Wilson loop W l = Ŵ l measured at the ground state is shown in Fig. 1 (d), where W y takes the exact quantized value at H = 0. After the time reversal symmetry breaking by the magnetic filed, W y takes nearly quantized values for a finite range of magnetic fields at H < H AFM * (H < H FM * ) for AFM (FM) Kitaev coupling. This indicates the Z 2 gauge structure in the pure S=1 Kitaev model remains a good description of perturbed Kitaev model away from the zero-field limit. Here we should note W y > 0 for L y = 3 and W y < 0 for L y = 4, while its absolute value is nearly quantized. This indicates distinct topological sectors can be accessed with changing L y . Beyond H FM * or H AFM * , W y tends to be negligibly small, particularly for wider systems, implying the absence of Z 2 gauge structure. The sudden drop in W y characterizes the phase transition at H AFM * and H FM * . The KSL Phase.-Below we show that the phase in Fig. 1(a) at 0 < H < H AFM * in AFM model or 0 < H < H FM * in FM model is a gapped KSL phase with topological nature. We first examine the spin-spin correlations S ij ≡ S i · S j − S i · S j and the Von Neumann entanglement entropy S vN to show the KSL is a gapped phase, and then we probe its topological properties from the Wilson loop,  the entanglement spectra and the edge modes. In the Kitaev limit (i.e., H=0), we confirm that the spin-spin correlations vanish beyond the nearest-neighbor links [52], which are similar to S = 1/2 pure Kitaev model [5].
where ρ A is the reduced density matrix of part A for the bipartition of the system into A and B, with ρ A obtained by tracing out the degrees of freedom of the B part. Here, we consider the cut parallel to e y and measure S vN for each cut with subsystem length L A . As shown in Fig. 2 (b), we find that S vN are independent on the positions of each cut and display flat behavior for cylinders with different widths, suggesting zero central charge or the existence of finite energy gap in the bulk.
In addition, we also examine the entanglement spectra for the cut parallel to e y at L A = L x /2, as shown in Fig. 2 (c), with a finite entanglement gap ∆ ξ separating the groundstate manifold from the higher spectra. On cylinders of width L y = 4, the entanglement spectra displays fourfold (quasi-)degeneracy [see Fig. 2 (c)], indicating the existence of boundary zero modes and its topological nature. The structure of entanglement spectra is similar to the topological KSL in spin S=1/2 Kitaev model [30]. The topological nature is also exhibited by the Wilson loop operatorŴ y , which takes nearly quantized mean value in the gapped KSL phase as shown in Fig. 1 (c). In addition, the distribution of W y along the cylinders is uniform with the nearly quantized value, as shown in the inset of Fig. 2 (d). All of these facts imply the Z 2 gauge structure remains a good description of the gapped KSL phase. To explicitly show the presence of edge modes, we compute the local on-site magnetization distributions S · e H [e H = (1, 1, 1)] along the zigzag chain of the cylinders, as shown in Fig. 2 (d), where the spins near the edge are much easier to be polarized than the bulk, indeed implying the gapless edge modes.
From the above analysis, the nearly quantized Wilson loop operator mean-value, the four-fold degeneracy in the entanglement spectra and the existence of edge modes all manifest the topological nature of the gapped KSL phase. Different from the FM Kitave model, the magnetic susceptibility χ in Fig. 1 (b) suggests the AFM Kitaev model further hosts an intermediate phase sandwiched between the gapped spin liquid and partially polarized paramagnet phase at H AFM * < H < H AFM * * . The intermediate phase is distinct from KSL phase by the vanishingly small Wilson loop, the disappearance of the degeneracy in entanglement spectra together with the edge mode. As shown in Fig. 3 (a) and its inset, S ij decays exponentially while S vN is flat. However, the local magnetization S · e H shows strong oscillation in real space, which is similar to the intermediate spin liquid phase in the S=1/2 Kitaev model. Due to the nonuniform pattern of magnetization, we suspect much larger system sizes are demanded to resolve the nature of the intermediate phase, which is beyond the scope of this work.
Instability of the KSL.-In the above, we have found that the gapped KSL in the AFM Kitaev model is much more robust against increasing magnetic field than the one in the FM Kitaev model (i.e., H AFM * H FM * ), though at zero field these two models are exactly equivalent [7]. For finite magnetic fields, the indications of such significant difference can be gained from the low-field magnitude of susceptibility χ, as shown in Fig. 1 (c). Our findings suggest a stable gapped KSL may be observable under applying suitable magnetic fields in materials harboring the AFM Kitaev interactions such as A 3 Ni 2 XO 6 (A=Li,Na, X=Bi,Sb) [58].
Furthermore, in the S=1 Kitaev materials, the superexchange processes would also give rise to the Heisenberg interactions between nearest-neighbor sites [58]. Thus we study the instability of the KSL against the Heisenberg interactions in the vicinity of the pure Kitaev model. Here we consider both AFM and FM Heisenberg coupling J. The phase transitions are characterized by the singular behavior of entanglement entropy and the kinks in the ground state energy density. As illustrated in the Figs. 4 (a-b) [56,58]. In other words, if the system is tuned proximate to the Kitaev phase at a finite J, the spin liquids we have discovered can be still induced with applying a magnetic field along [111] direction.
Discussions and Summary.-In this work, we investigate the magnetic field induced phases in the S=1 Kitaev honeycomb model by extensive large scale DMRG simulations. In the weak field limit, we identify a gapped KSL phase with topological nature for both the AFM and FM Kitaev couplings, including Z 2 gauge structure, the edge modes and fourfold degenerate entanglement spectra on cylinders. Interestingly, such gapped KSL phase is much more stable against increasing magnetic field for the AFM Kitaev coupling than the FM Kitaev coupling, which resembles to the spin S=1/2 Kitaev model subject to a perpendicular magnetic field [30]. Furthermore, when the gapped KSL is destroyed by the large field in the AFM Kitaev model, we also find an intermediate nonmagnetic phase before entering the partial polarized phase, which may be a critical phase similar to spin S=1/2 case. Theoretical understanding of spin S=1 Kitaev model in a magnetic field is still a challenge partially due to the lack of exact solution even in the pure Kitaev limit, and thus it will be an excellent task to pursue in the future. The recently proposed S=1 Kitaev materials such as layered transition metal oxides A 3 Ni 2 XO 6 (A=Li,Na, X=Bi,Sb) [58] host AFM Kitaev coupling with honeycomb plane perpendicular to [111] direction, thus that our findings can be ideally tested in these systems. Our work may also inspire the search for the spin liquids in possible S=1 and high-S Kitaev materials in the fu- ture. The signatures of the spin liquids in the spin-1 Kitaev model such as the quantized thermal Hall conductivity and fractionalized excitations in inelastic neutron scattering can be directly probed, which is technically similar to the intensively investigated S=1/2 Kitaev materials. Note added: At the final stage of finishing this work, we noticed that Ref. [57] constructed a tensor network groundstate wavefunction for the pure S=1 Kitaev model. The magnetic field effect was also examined in Ref. [57] via the magnetizations, which are consistent with our results in Fig. 1  (b-c). In our work, we identify the magnetic field induced AFM/FM KSL from different aspects of systematical characterization based on unbiased DMRG method. We should also note that the largest field for the AFM model in Ref. [57] is up to H ≈ 0.3/ √ 3 < H AFM * , which lies within the AFM KSL phase identified in our work.
Acknowledgments-We acknowledge the computational resources at CSUN and KITS for performing the numerical simulations in this work. This material is