Anyon polarons as a window into the competing phases of the Kitaev-Gamma-Gamma' model

TL;DR

This paper explores how anyon gap closures in the Kitaev honeycomb model with perturbations reveal phase transitions into various magnetic orders and spin-liquid states, providing insights into the model's complex phase diagram.

Contribution

It introduces a method to understand phase transitions via anyon gap closures in the Kitaev-Gamma-Gamma' model, linking them to specific competing phases and critical coupling values.

Findings

Anyon gap closings signal phase transitions into zigzag, stripy, 120°, and incommensurate spiral phases.

In certain perturbations, the spin liquid coexists with magnetic order due to selective gap closures.

The study aligns theoretical gap-closing points with numerical phase boundaries.

Abstract

We investigate the dispersions of anyon quasi-particles in the Kitaev honeycomb spin-liquid perturbed by $\Gamma$ and $\Gamma'$ couplings in order to understand phase transitions into competing states through anyon gap-closing instabilities. We demonstrate how anyon gap closings allow to understand phase transitions into a plethora of previously identified competing phases -- including zigzag, stripy, $120^\circ$, and incommensurate spiral phases -- and are in agreement with numerical studies not only on the nature of the phases, but also on the specific critical values of $\Gamma$ and $\Gamma'$ couplings. Remarkably, when the anti-ferromagnetic Kitaev model is perturbed by an ferromagnetic $\Gamma$ interaction, we find that the single-vison and fermion gaps remain open while the gap of a magnon-like local boson vanishes, implying that the resulting state has coexistence of a spontaneous broken symmetry and the fractionalization pattern of the Kitaev spin liquid. The magnetic long-range order could be either a stripy antiferromagnet or an incommensurate spiral, depending on the sign of $\Gamma'$.

Equations

$$H_{⟨i,j{⟩}_{\alpha }}=$$

$$+{\Gamma }^{\prime }({\sigma }_i^{\alpha }{\sigma }_j^{\beta }+{\sigma }_i^{\alpha }{\sigma }_j^{\gamma }+{\sigma }_i^{\beta }{\sigma }_j^{\alpha }+{\sigma }_i^{\gamma }{\sigma }_j^{\alpha }).$$

$$H_{\mathrm{v}}=\sum _{⟨i,j⟩}{t}_{i,j}^{{\Gamma }^{\prime }}|v_i⟩⟨v_j|+\sum _{⟨⟨i,j⟩⟩}{t}_{i,j}^{\Gamma }|v_i⟩⟨v_j|+\sum _i{\Delta }_v|v_i⟩⟨v_i|.$$

$${\tilde{\chi }}_{\mathbf{r},\alpha }^{†}|\Omega ⟩={\chi }_{\mathbf{r},\alpha }^{†}|0;{\Psi }_c(\mathbf{r},\alpha )⟩.$$

$$H_{\mathrm{f}}=H_c(0)+\sum _{\mathbf{r}\in A,\alpha }{\Delta }_{\tilde{\chi }}{\tilde{\chi }}_{\mathbf{r},\alpha }^{†}{\tilde{\chi }}_{\mathbf{r},\alpha }+H_{\mathrm{f}}^{\Gamma }+H_{\mathrm{f}}^{{\Gamma }^{\prime }}.$$

$$H_{\mathrm{f}}^{\Gamma }$$

$$+T_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{\Gamma ,yx}{\tilde{\chi }}_{\mathbf{r},y}^{†}{\tilde{\chi }}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,x}+T_{{\mathbf{a}}_1}^{\Gamma ,zy}{\tilde{\chi }}_{\mathbf{r},z}^{†}{\tilde{\chi }}_{\mathbf{r}-{\mathbf{a}}_1,y}$$

$$+T_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{\Gamma ,xz}{\tilde{\chi }}_{\mathbf{r},x}^{†}{\tilde{\chi }}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,z}+T_{-{\mathbf{a}}_2}^{\Gamma ,yx}{\tilde{\chi }}_{\mathbf{r},y}^{†}{\tilde{\chi }}_{\mathbf{r}+{\mathbf{a}}_2,x}+H.c.],$$

$$H_{\mathrm{f}}^{{\Gamma }^{\prime }}$$

$$+T_{-{\mathbf{a}}_2}^{{\Gamma }^{\prime },zx}{\tilde{\chi }}_{\mathbf{r},z}^{†}{\tilde{\chi }}_{\mathbf{r}+{\mathbf{a}}_2,x}+T_0^{{\Gamma }^{\prime },yz}{\tilde{\chi }}_{\mathbf{r},y}^{†}{\tilde{\chi }}_{\mathbf{r},z}$$

$$+T_0^{{\Gamma }^{\prime },xy}{\tilde{\chi }}_{\mathbf{r},x}^{†}{\tilde{\chi }}_{\mathbf{r},y}+T_0^{{\Gamma }^{\prime },zx}{\tilde{\chi }}_{\mathbf{r},z}^{†}{\tilde{\chi }}_{\mathbf{r},x}+H.c.]$$

$$+\sum _{\mathbf{r}\in A,\alpha }\sum _{\mathbf{R}}{p}_{\mathbf{R},\alpha }(c_{\mathbf{r}+\mathbf{R}}-i{c}_{\mathbf{r}+{\mathbf{\delta }}_{\alpha }-\mathbf{R}}){\tilde{\chi }}_{\mathbf{r},\alpha }+H.c.$$

$$T_{-{\mathbf{a}}_2}^{\Gamma ,zy}$$

$$=-\Gamma ⟨0;{\Psi }_c(\mathbf{r},z)|$$

$$\left(1+i{c}_{\mathbf{r}+{\mathbf{a}}_2}{c}_{\mathbf{r}+{\mathbf{\delta }}_z}\right)|0;{\Psi }_c(\mathbf{r}+{\mathbf{a}}_2,y)⟩,$$

$$T_0^{{\Gamma }^{\prime },yz}$$

$$+H_{\mathbf{r}-{\mathbf{a}}_1,y}^{{\Gamma }^{\prime }}+H_{\mathbf{r}-{\mathbf{a}}_2,z}^{{\Gamma }^{\prime }}]{\tilde{\chi }}_{\mathbf{r},z}^{†}|\Omega ⟩$$

$$={\Gamma }^{\prime }⟨0;{\Psi }_c(\mathbf{r},y)|c_{\mathbf{r}}g|0;{\Psi }_c(\mathbf{r},z)⟩,$$

$$g$$

$$T_{-{\mathbf{a}}_2}^{\Gamma ,zy}$$

$$=T^{\Gamma }\in \mathbb{R},$$

$$T_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{{\Gamma }^{\prime },yz}$$

$$=T^{{\Gamma }^{\prime }}\in \mathbb{R},$$

$$d_{\mathbf{r},\alpha ,l}^{†}|\Omega ⟩={\chi }_{\mathbf{r},\alpha }^{†}{\alpha }_{(\mathbf{r},\alpha ),l}^{†}|0;{\Psi }_c(\mathbf{r},\alpha )⟩.$$

$$o_{\mathbf{r},\alpha }^{†}|\Omega ⟩={\chi }_{\mathbf{r}+{\mathbf{\delta }}_{\alpha }-{\mathbf{\delta }}_{\gamma },\gamma }^{†}{\chi }_{\mathbf{r},\beta }^{†}|0;{\Psi }_c(\mathbf{r}+{\mathbf{\delta }}_{\alpha }-{\mathbf{\delta }}_{\gamma },\gamma ;\mathbf{r},\beta )⟩.$$

$$H_{\mathrm{b}}^d$$

$$+\sum _{\mathbf{r}\in A}[M_{-{\mathbf{a}}_2}^{\Gamma ,zy}{d}_{\mathbf{r},z}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_2,y}+M_{{\mathbf{a}}_1}^{\Gamma ,xz}{d}_{\mathbf{r},x}^{†}{d}_{\mathbf{r}-{\mathbf{a}}_1,z}$$

$$+M_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{\Gamma ,yx}{d}_{\mathbf{r},y}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,x}+M_{{\mathbf{a}}_1}^{\Gamma ,zy}{d}_{\mathbf{r},z}^{†}{d}_{\mathbf{r}-{\mathbf{a}}_1,y}$$

$$+M_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{\Gamma ,xz}{d}_{\mathbf{r},x}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,z}+M_{-{\mathbf{a}}_2}^{\Gamma ,yx}{d}_{\mathbf{r},y}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_2,x}$$

$$+M_{-{\mathbf{a}}_1+{\mathbf{a}}_2}^{{\Gamma }^{\prime },yz}{d}_{\mathbf{r},y}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,z}+M_{{\mathbf{a}}_1}^{{\Gamma }^{\prime },xy}{d}_{\mathbf{r},x}^{†}{d}_{\mathbf{r}-{\mathbf{a}}_1,y}$$

$$+M_{-{\mathbf{a}}_2}^{{\Gamma }^{\prime },zx}{d}_{\mathbf{r},z}^{†}{d}_{\mathbf{r}+{\mathbf{a}}_2,x}+M_0^{{\Gamma }^{\prime },yz}{d}_{\mathbf{r},y}^{†}{d}_{\mathbf{r},z}$$

$$+M_0^{{\Gamma }^{\prime },xy}{d}_{\mathbf{r},x}^{†}{d}_{\mathbf{r},y}+M_0^{{\Gamma }^{\prime },zx}{d}_{\mathbf{r},z}^{†}{d}_{\mathbf{r},x}+H.c.],$$

$$H_{\mathrm{b}}^o$$

$$+\sum _{\mathbf{r}\in A}[t_{-{\mathbf{a}}_1}^{o,zx}{o}_{\mathbf{r},z}^{†}{o}_{\mathbf{r}+{\mathbf{a}}_1,x}+t_{-{\mathbf{a}}_1}^{o,yz}{o}_{\mathbf{r},y}^{†}{o}_{\mathbf{r}+{\mathbf{a}}_1,z}$$

$$+t_{{\mathbf{a}}_1-{\mathbf{a}}_2}^{o,xy}{o}_{\mathbf{r},x}^{†}{o}_{\mathbf{r}-{\mathbf{a}}_1+{\mathbf{a}}_2,y}+t_{{\mathbf{a}}_1-{\mathbf{a}}_2}^{o,zx}{o}_{\mathbf{r},z}^{†}{o}_{\mathbf{r}-{\mathbf{a}}_1+{\mathbf{a}}_2,x}$$

$$+t_{{\mathbf{a}}_2}^{o,yz}{o}_{\mathbf{r},y}^{†}{o}_{\mathbf{r}-{\mathbf{a}}_2,z}+t_{{\mathbf{a}}_2}^{o,xy}{o}_{\mathbf{r},x}^{†}{o}_{\mathbf{r}-{\mathbf{a}}_2,y}+H.c.],$$

$$H_{\mathrm{b}}^{do}$$

$$+{\mathcal{T}}_4^x{d}_{\mathbf{r}-{\mathbf{a}}_2,z})+o_{\mathbf{r},y}^{†}({\mathcal{T}}_1^y{d}_{\mathbf{r},z}+{\mathcal{T}}_2^y{d}_{\mathbf{r},x}$$

$$+{\mathcal{T}}_3^y{d}_{\mathbf{r}+{\mathbf{a}}_1-{\mathbf{a}}_2,z}+{\mathcal{T}}_4^y{d}_{\mathbf{r}+{\mathbf{a}}_1,x})+o_{\mathbf{r},z}^{†}({\mathcal{T}}_1^z{d}_{\mathbf{r},x}$$

$$+{\mathcal{T}}_2^z{d}_{\mathbf{r},y}+{\mathcal{T}}_3^z{d}_{\mathbf{r}+{\mathbf{a}}_2,x}+{\mathcal{T}}_4^z{d}_{\mathbf{r}-{\mathbf{a}}_1+{\mathbf{a}}_2,y})+H.c.$$

$$d_{\mathbf{r},\alpha }^{†}=(d_{\mathbf{r},\alpha ,1}^{†},\dots ,d_{\mathbf{r},\alpha ,N}^{†}).$$

$$[M_{-{\mathbf{a}}_2}^{\Gamma ,zy}{]}_{m,n}$$

$$=-\Gamma ⟨0;{\Psi }_c(\mathbf{r},z)|{\alpha }_{(\mathbf{r},z),m}\left[1+i{c}_{\mathbf{r}+{\mathbf{a}}_2}{c}_{\mathbf{r}+{\mathbf{\delta }}_z}\right]$$

$${\alpha }_{(\mathbf{r}+{\mathbf{a}}_2,y),n}^{†}|0;{\Psi }_c(\mathbf{r}+{\mathbf{a}}_2,y)⟩,$$

$$[M_0^{{\Gamma }^{\prime },yz}{]}_{m,n}$$

Full text

I Introduction

Quantum spin liquids (QSLs) are highly entangled states of matter which feature the emergence of fractionalized excitations [1, 2, 3, 4]. Among the various theoretical models proposed to host QSLs [5, 6, 7, 8, 9], the Kitaev honeycomb model [9] has emerged as a paradigmatic example due to its exact solvability and potential realizations [10] in Kitaev materials such as Na2IrO3, $\alpha$-Li2IrO3, and $\alpha$-RuCl3 [11, 12, 13]. While the ideal Kitaev model realizes a $\mathbb{Z}_{2}$ QSL with itinerant Majorana fermion (spinon) and static $\mathbb{Z}_{2}$ gauge flux (vison) excitations, non-Kitaev interactions present in Kitaev materials—such as the Heisenerg exchange ($J$) [10], off-diagonal $\Gamma$ and $\Gamma^{\prime}$ terms [14, 15, 16], and external magnetic fields—can destabilize the Kitaev spin liquid (KSL) and drive the system into other states [14, 17, 18, 19, 20, 21]. In particular, $\alpha$-RuCl3 is believed to be described by a nearest-neighbor $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model, possibly supplemented by third-neighbor Heisenberg interactions ($J_{3}$) and inter-layer couplings, although accurately determining these parameters remains a challenge [22, 23, 24, 25]. These additional couplings are believed to be responsible for stabilizing the observed zigzag magnetic order at low temperatures, which itself can be suppressed by external magnetic fields and replaced by a novel phase—whose nature remains under active debate [26, 27, 28, 29, 30, 31]. Therefore, understanding the impact of these additional couplings on the ideal KSL remains a major open problem, the resolution of which could help elucidate more clearly the connections between extended Kitaev models and real materials.

Numerical studies of such extended Kitaev models have uncovered rich phase diagrams featuring both magnetically ordered phases and other competing novel QSLs (see Refs. [11, 32] for reviews of the results). When analyzing extended Kitaev models harboring multiple competing orders, a powerful starting point is to regard the more symmetric KSL as the mother state, with broken-symmetry phases—such as magnetically ordered states—emerging as its descendants (for earlier studies exploiting this philosophy, see, e.g., Refs. [33, 34]).

A large class of transitions from the parent QSL to nearby phases can be understood as instabilities where the gap for some of its quasiparticles closes. In particular, if a QSL hosts a quasiparticle with bosonic self-statistics, its gap closing (i.e., condensation) gives rise to a new phase in which any other quasiparticle with nontrivial mutual statistics becomes “confined”, namely it no longer exists as a finite-energy excitation (see, e.g., Ref. [35]). For instance, in the toric code $\mathbb{Z}_{2}$ QSL, the condensation of one of its non-local bosons (e.g. the $e$-particle) confines the other anyons, resulting in a conventional magnetically ordered state (see, e.g., Ref. [36]). By contrast, if the condensed quasiparticle is a local boson (i.e., an ordinary magnon-like quasiparticle), which has trivial statistics with all quasiparticles but carries non-trivial symmetry quantum numbers, its gap closing results in a new phase with the same anyon content of the parent QSL, but with spontaneously broken symmetry dictated by its quantum numbers. In the case of the ideal Kitaev model, there exists three quasiparticle anyon types [9]: local bosons, fermions, and a flux particle also known as vison. Each vison can be viewed as a bare boson dressed by a cloud of fermions since it acts as a $\pi$-flux for the latter. This “cloud” is responsible for transforming the bare bosonic statics of the vison and rendering it, for instance, a non-Abelian anyon when the fermions have a Bogoliubov Chern band in the presence of a Haldane mass term [9, 37, 38]. The gap closing of the single visons is an indication of confinement, which can be viewed as resulting from the condensation of their bare bosonic core 111We are thankful to Chao-Ming Jian for helping us clarify this point. Moreover, a pair of visons can form a bound state which can be either a local boson or a fermion. Therefore, if the lowest-energy quasiparticle is a bosonic vison-pair bound-state and its gap closes under a perturbation without the accompanying gap closing of single visons, one expects a phase transition from the KSL into a symmetry-broken state that nonetheless retains the same anyon content as the KSL. In this case the pattern of symmetry breaking can be determined from the quantum of the condensed bosonic vison-pair.

As will be detailed in latter sections, by computing the dispersions and the quantum numbers of bosonic vison pairs of KSL, we will be able to determine a variety of magnetic orders observed in extended Kitaev models, including zigzag, stripy, $120^{\circ}$ (FM${}_{U_{6}}$ and AFM${}_{U_{6}}$, see definition below) [40, 41]. In most of these cases, we observe that the single visons are also condensed at nearby values of the additional couplings, indicating that many of the resulting phases are likely trivial magnetically ordered states. However, in the anti-ferromagnetic Kitaev model with negative $\Gamma$, we find a regime where the bosonic visons-pairs condense while single visons remain gapped. The resulting phase is therefore expected to be a QSL with the same anyon content of the KSL but coexisting with long-range magnetic order, which we find to be either stripy or incommensurate (IC) order depending on the sign of $\Gamma^{\prime}$. Notably, an earlier study by Zhang et al. [40] employed a similar vison-pair approach to analyze the $K$–$J$–$\Gamma$ model and found results consistent with our findings regarding the $\Gamma$-induced orders, such as the FM${}_{U_{6}}$ ($\Gamma>0$) and IC ($\Gamma<0$) orders in the anti-ferromagnetic Kitaev model. However, the single-vison behavior was not considered in this study, leaving open the question of whether the resulting phases are trivial magnetically ordered states or QSLs with coexisting symmetry breaking, as we are finding. We hope that future numerical studies could identify these interesting phases that combine symmetry breaking with fractionalization, by investigation of both local magnetic orders and entanglement signatures [42, 43], or other indicators such as topological degeneracy in the torus.

Moreover, since bosonic vison pairs are created by local operators, they play a central role in understanding typical experimentally measurable correlation functions of local operators, such as the dynamical spin structure factor $S(\omega,\bm{Q})$ [40]. Owing to their bosonic character, the bosonic vison pairs can naturally capture the emergence of magnon-like low-energy dispersions in $S(\omega,\bm{Q})$ near the development of magnetic long-range orders. Finally, we note that the fermionic spinons of the KSL may also undergo gap-closing/opening transitions accompanied by changes in their band Chern numbers, which in turn can modify the nature of vison excitations [9, 41, 44]. We therefore also analyze in detail the gap closings associated with these fermionic quasiparticles.

I.1 The model and some previous results

Analyzing the dispersion of anyons has provided valuable insights into extended Kitaev models, such as the Kitaev model with a Zeeman field [38, 41, 44] and the $K$–$J$–$\Gamma$ model [40]. In this work, we further apply it to the $K$–$\Gamma$–$\Gamma\textquoteright$ model [20, 45, 21, 46] to assess its consistency with numerics and to extract new physical insights. The model is given by $H_{K\text{-}\Gamma\text{-}\Gamma^{\prime}}=\sum_{\langle i,j\rangle_{\alpha}}H_{\langle i,j\rangle_{\alpha}}$, with

[TABLE]

Here $\langle i,j\rangle_{\alpha}$ denotes a bond of type $\alpha=x,y,z$, and $\beta\neq\gamma\neq\alpha$. Previous numerical studies have revealed several important features of this model. In the ferromagnetic (FM) Kitaev model [$K=-1$ in Eq. (I.1)], which is considered pertinent to $\alpha$-RuCl3, the KSL is highly fragile against the $\Gamma$ interaction. When $\Gamma\gtrsim 0.03$, the KSL has been argued to be replaced by a novel QSL whose nature remains elusive: tensor product state and density matrix renormalization group (DMRG) calculations suggest a nematic paramagnetic (NP) phase characterized by anisotropic bond energies [21, 47], whereas variational Monte Carlo calculations find a proximate KSL phase featuring multiple Majorana cones [48, 49]. Both phases, however, are destabilized in favor of zigzag order upon the inclusion of a small negative $\Gamma^{\prime}\approx-0.03$ [21, 49]. For $\Gamma<0$, the KSL is replaced by the AFM${}_{U_{6}}$ order (referred to as RS${}_{U_{6}}$ in Ref. [46]), which maps to a Néel order under a six-site unitary transformation $U_{6}$ [50, 51, 46, 32]. On the other hand, the KSL is more robust against the $\Gamma$ interaction in the antiferromagnetic (AFM) Kitaev model [$K=1$ in Eq. (I.1)]. When $\Gamma\gtrsim 0.3$, the KSL is replaced by a FM${}_{U_{6}}$ order which maps to a FM order under the $U_{6}$ transformation [46, 51, 32]. In contrast, for $\Gamma\lesssim-0.5$, the system could host either stripy or IC spiral order [14, 46]. The role of the $\Gamma^{\prime}$ interaction in this regime is less studied. The goal of this study is to provide a microscopic understanding of how various phases emerge from the KSL in the $K$-$\Gamma$-$\Gamma^{\prime}$ model by starting from the ideal Kitaev limit ($\Gamma=\Gamma^{\prime}=0$), computing perturbatively in $\Gamma$ and $\Gamma^{\prime}$ the dispersions of its anyon quasiparticles, and investigating the resulting anyon gap-closing instabilities.

II Quasiparticles of the KSL

Owing to its exact solvability, the Kitaev honeycomb model can be adopted as a well-controlled starting point. One can therefore start from the KSL and explore, in a perturbative manner, the impact of the $\Gamma$ and $\Gamma\textquoteright$ interactions on its quasiparticle excitations [40, 44, 52, 38, 41]. In this work, we focus primarily on the leading (linear) order effects of the $\Gamma$ and $\Gamma^{\prime}$ interactions. The KSL contains three types of quasiparticles: single visons, fermions including itinerant Majorana fermions and fermionic vison pairs, and local bosonic vison pairs. An important effect of the non-Kitaev interactions like $\Gamma$ and $\Gamma^{\prime}$ is that they induce fluctuations in the gauge degrees of freedom, rendering single visons and vison pairs—both fermionic and bosonic—dynamical, and thereby making the model no longer exactly solvable. In the following, we analyze the impact of the $\Gamma$ and $\Gamma^{\prime}$ interactions on each type of quasiparticle in the KSL.

II.1 Single visons

The single visons move on an effective triangular lattice formed by the hexagons of the honeycomb lattice (see Fig. 1(c)) [52, 38]. The $\Gamma$/$\Gamma^{\prime}$ interaction enables a single vison to hop between next-nearest-neighbor (NNN)/nearest-neighbor (NN) hexagons. To linear order in $\Gamma$/$\Gamma^{\prime}$, due to the existence of a vison gap $\Delta_{v}$, the effective Hamiltonian of a single vison only contains hoppings within the degenerate states associated with its position 222Since physical states only contain pairs of visons, one imagines a distant auxiliary second vison which is not affected by the perturbation.. Therefore, the effective Hamiltonian for single visons has the form:

[TABLE]

Here $|v_{i}\rangle$ is the state with a vison at plaquette $i$. $t_{i,j}^{\Gamma/\Gamma^{\prime}}$ is the $\Gamma$/$\Gamma^{\prime}$-induced hopping amplitude, defined as the matrix element of the $\Gamma/\Gamma^{\prime}$ interaction between states where a vison locates at sites $i$ and $j$ (see its expression in Supplemental Material). $\Delta_{v}\approx 0.15|K|$ is the single-vison creation energy [9].

Interestingly, following the computation methods previously developed in Refs. [37, 38, 41, 52], we have found that for the FM Kitaev model, $t^{\Gamma}/\Gamma$ is finite while $t^{\Gamma^{\prime}}/\Gamma^{\prime}=0$; For the AFM Kitaev model, we find that $t^{\Gamma}/\Gamma=0$ whereas $t^{\Gamma^{\prime}}/\Gamma^{\prime}$ is finite, in agreement with Ref. [52] (see Fig. 2). Moreover, our computation reveals that the $t^{\Gamma^{\prime}}_{i,j}$ values give rise to a non-trivial projective implementation of translational symmetry, corresponding to a $\pi$-flux per honeycomb unit cell experienced by visons in the AFM Kitaev model, thereby effectively doubling their unit cell (see vison bands in Fig. 3). As we have shown previously [37, 38, 41], this pattern of translational symmetry fractionalization is an intrinsic property of the AFM-KSL, independent of the perturbation that induces vison mobility, making it a sharply distinct symmetry enriched topologically ordered state compared to the FM-KSL, whose visons have a trivial implementation of lattice translations. The significant (vanishing) value of $t^{\Gamma}/\Gamma$ in the FM (AFM) Kitaev model implies that the single-vison gap closes at relatively small (large) $\Gamma$, consistent with the weaker (stronger) stability of KSL against $\Gamma$ interaction in the FM (AFM) case [14, 40]. We have also examined the second-order single-vison hopping $t^{\Gamma,(2)}$ induced by $\Gamma$, and found a similar qualitative behavior: $t^{\Gamma,(2)}/\Gamma^{2}\neq 0$ ($=0$) for FM (AFM) Kitaev model (see Supplemental Material for details). This supports the idea that, in the AFM Kitaev model, single visons remain gapped—that is, matter fermions remain deconfined—until much larger $\Gamma$ values, opening the door to a novel $\Gamma$-induced QSL phase exhibiting spontaneous symmetry breaking (see discussions in later sections). Although the linear-order single-vison hopping induced by the $\Gamma^{\prime}$ interaction vanishes in the FM Kitaev model, higher-order contributions are expected to render it finite.

Crutially, unlike the case of a Zeeman field—which induces a Haldane mass gap on the matter fermions [7], and renders the vison an exponentially localized quasiparticle ($\sigma$ anyon of the Ising topological order)—the $\Gamma$ and $\Gamma^{\prime}$ interactions do not immediately gap out the fermions. In the present context, this implies that visons behaves as vortices which are dressed by power-law decaying modification of the surrounding matter fermion fluid. Therefore, we need to exercise caution to correctly extract their hopping amplitude numerically since the states in a torus always contain at least two visons. To do so, we have introduced an artificial Haldane mass ($\kappa$) term to gap out the matter fermions, allowing us to define visons as well-separated individual excitations. We then extract the thermodynamic-limit values of $t^{\Gamma}/\Gamma$ and $t^{\Gamma^{\prime}}/\Gamma^{\prime}$ by taking the limit $\kappa\rightarrow 0$, while keeping the vison–vison separation remains greater than their effective localization length $\xi_{c}=v_{c}/\Delta_{c}$, where $v_{c}=\sqrt{3}|K|$ is the $c$-fermion Fermi velocity and $\Delta_{c}=6\sqrt{3}\kappa$ is the gap induced by the mass term.

In contrast to single visons, vison pairs (including two NN or NNN visons) induced a power-law modification of the surrounding matter fermion fluid that decays with a larger power with distance [40], leading to a faster convergences for the estimate of their hopping amplitudes with system size. Below, we will analyze their properties in the presence of the $\Gamma$ and $\Gamma^{\prime}$ interactions.

II.2 Fermions

A fermionic vison pair consists of two NN visons (see the illustration in Fig. 1(c)) accompanied by a modified matter fermion background [54, 55, 40, 44, 41]. Utilizing the Majorana representation of local spins ($\sigma_{\bm{r}}^{\alpha}\leftrightarrow ib_{\bm{r}}^{\alpha}c_{\bm{r}}$) [9], a fermionic vison pair on bond $(\bm{r},\alpha)$ with $\bm{r}\in A$ is defined as 333The projection operator $P=\prod_{\bm{r}}(1+D_{\bm{r}})/2$ and the normalization factor $\mathcal{N}=2^{(N_{s}-1)/2}$ have been omitted, same for the bosonic vison pairs. $D_{\bm{r}}=b_{\bm{r}}^{x}b_{\bm{r}}^{y}b_{\bm{r}}^{z}c_{\bm{r}}$ and $N_{s}$ is the number of sites in the system. Here we are considering a periodic system (torus).:

[TABLE]

Here $|\Omega\rangle=|0;\Psi_{c}(0)\rangle$ denotes the ground state of the Kitaev Hamiltonian without any flux. $\chi_{\bm{r},\alpha}^{\dagger}=b_{\bm{r}}^{\alpha}-ib_{\bm{r}+\bm{\delta}_{\alpha}}^{\alpha}$ is the creation operator of a $\mathbb{Z}_{2}$ gauge fermion defined on bond $(\bm{r},\alpha)$. The state $|0,\Psi_{c}(\bm{r},\alpha)\rangle$ is the ground state of a matter $c$-fermion BdG Hamiltonian $H_{c}(\bm{r},\alpha)$ incorporating a pair of $\pi$ fluxes separated by the bond $(\bm{r},\alpha)$. Importantly, $|0,\Psi_{c}(\bm{r},\alpha)\rangle$ has the same (even) total fermion parity as $|\Omega\rangle$, so $\tilde{\chi}_{\bm{r},\alpha}^{\dagger}|\Omega\rangle$ carries the opposite (odd) fermion parity and is therefore referred to as a fermionic vison pair [57, 55, 40, 38].

The $\Gamma$ interaction induces hopping of $\tilde{\chi}$-fermions [40], while the $\Gamma^{\prime}$ interaction not only generates additional $\tilde{\chi}$-hopping but also leads to hybridization between the $\tilde{\chi}$ and $c$ fermions (see Fig. 4(c)). Specifically, the $\Gamma^{\prime}$ interaction gives rise to tunneling as well as pair creation and annihilation of $\tilde{\chi}$ and $c$ fermions (see Supplemental Material). The effective Hamiltonian of the fermion quasiparticles is therefore:

[TABLE]

Here $H_{c}(0)$ is the $c$-fermion BdG Hamiltonian in the zero-flux sector, $\Delta_{\tilde{\chi}}\approx 0.26|K|$ is the excitation gap of a fermion vison pair in Kitaev model [41, 40]. The part induced by $\Gamma$ and $\Gamma^{\prime}$ interactions reads:

[TABLE]

The fermion hopping parameter $T^{\Gamma/\Gamma^{\prime},\alpha\beta}_{-\bm{R}}$ ($\bm{R}$ is a lattice vector) is defined as the matrix element of the $\Gamma/\Gamma^{\prime}$ interaction between states $\langle\Omega|\tilde{\chi}_{\bm{r},\alpha}$ and $\tilde{\chi}_{\bm{r}+\bm{R},\beta}^{\dagger}|\Omega\rangle$. E.g.,

[TABLE]

[TABLE]

Here $H^{\Gamma/\Gamma^{\prime}}_{\bm{r},\alpha}$ stands for the $\Gamma/\Gamma^{\prime}$ interaction on bond $(\bm{r},\alpha)$. Note for each $H^{\Gamma^{\prime}}_{\bm{r},\alpha}$ in Eq. (II.2), there is only one two-spin term contributes to the hopping process. Due to rotational symmetry, the $\Gamma$-induced $\tilde{\chi}$-fermion hopping amplitudes satisfy:

[TABLE]

and $T^{\Gamma}/\Gamma\approx-1.44\ (-0.17)$ in the FM (AFM) Kitaev model (see Ref. [40] for plots of $T^{\Gamma}/\Gamma$ at different system sizes). Similarly, these same symmetries constraint the $\Gamma^{\prime}$ terms as follows:

[TABLE]

with $T^{\Gamma^{\prime}}/\Gamma^{\prime}\approx-0.27\ (0.61)$ in the FM (AFM) Kitaev model (see Figs. 4(a)–(b)). Note that the $\Gamma$-induced fermion hopping in the FM Kitaev model is an order of magnitude stronger than in the AFM Kitev model. This suggests that the fermions are more likely to undergo gap closing at much smaller $\Gamma$ values in the FM Kitaev model, consistent with the greater fragility of the KSL under this perturbation. By contrast, the $\Gamma^{\prime}$-induced fermion hopping is approximately twice as strong in the AFM model compared to the FM one.

While $T^{\Gamma}$ and $T^{\Gamma^{\prime}}$ exhibit strong dependence on the sign of $K$, the $\tilde{\chi}$-$c$ coupling $p_{\bm{R},\alpha}$ is found to be independent of it. Moreover, a $\tilde{\chi}_{\bm{r},\alpha}$ fermion couples only to nearby $c$ fermions, with $\max(|p_{\bm{R},\alpha}/\Gamma^{\prime}|)\approx 0.88$ (see Supplemental Material).

II.3 Bosons

A NN bosonic vison pair at bond $(\bm{r},\alpha)$ consists of a fermion vison pair and a matter fermion excitation [40, 38]:

[TABLE]

Here $\alpha_{(\bm{r},\alpha),l}^{\dagger}$ is the creation operator for the $l$-th ($l=1,\dots,N$, with $N$ being the numeber of unit cells in the system) Bogoliubov quasiparticle of the matter fermion Hamiltonian $H_{c}(\bm{r},\alpha)$, and $d_{\bm{r},\alpha,l}^{\dagger}$ denotes the creation operator for the $l$-th NN bosonic vison pair at bond $(\bm{r},\alpha)$. An important feature of the vison pair $d_{\bm{r},\alpha,l}^{\dagger}$ is that it can be created from, or annihilated into, the KSL vacuum $|\Omega\rangle$ by local operators, such as Pauli operators $\sigma_{\bm{r}}^{\alpha}$ and $\sigma_{\bm{r}+\bm{\delta}_{\alpha}}^{\alpha}$ or their products. These quasiparticles can be therefore viewed as generalized magnons whose condensation can induce spontaneous symmetry breaking and magnetic orders, as formulated in Sec. VI.1 [38, 40]. They also play an important role in the behavior of correlation functions of local operators, such as the dynamical spin structure factor $S(\omega,\bm{Q})$ (see its expression in Sec. VI.2) [54, 55, 40]. As $S(\omega,\bm{Q})$ in the ideal Kitaev model is quantitatively well reproduced by the bosonic vison-pair approximation (see plots in Fig. 10) [40, 55], extending this framework to extended Kitaev models can provide valuable insights into emergent magnetic orders driven by non-Kitaev interactions.

In addition to the NN bosons, we also consider the NNN bosonic vison pairs, as they have slightly lower energy than the NN vison pairs in the ideal Kitaev model (with excitation energy $\Delta_{o}\approx 0.24|K|$) [9]. A NNN vison pair with two visons being connected by the bond $(\bm{r},\alpha)$ is defined as (see illustration in Fig. 1(c)):

[TABLE]

Here $o_{\bm{r},\alpha}^{\dagger}$ is the creation operator for the NNN boson vison pair at bond $(\bm{r},\alpha)$. In Eq. (11), the indices $\beta$ and $\gamma$ are fixed by $\alpha$ through cyclic ordering. That is, $(\alpha,\beta,\gamma)$ always form an ordered cyclic permutation of $(x,y,z)$; e.g., if $\alpha=z$, then $\beta=x$ and $\gamma=y$. The state $|0;\Psi_{c}(\bm{r}+\bm{\delta}_{\alpha}-\bm{\delta}_{\gamma},\gamma;\bm{r},\beta)\rangle$ is the ground state of matter fermion Hamiltonian $H_{c}(\bm{r}+\bm{\delta}_{\alpha}-\bm{\delta}_{\gamma},\gamma;\bm{r},\beta)$, which hosts two $\pi$ fluxes connected by $(\bm{r},\alpha)$. $o_{\bm{r},\alpha}^{\dagger}|\Omega\rangle$ is bosonic since it shares the same fermion parity as $|\Omega\rangle$.

The $\Gamma$ interaction induces hopping of both NN and NNN bosons, whereas the $\Gamma^{\prime}$ interaction generates hopping of NN bosons and tunneling between NN and NNN bosons, as illustrated in Figs. 4(c)–(d). The boson effective Hamiltonian reads: $H_{\mathrm{b}}=H_{\mathrm{b}}^{d}+H_{\mathrm{b}}^{o}+H_{\mathrm{b}}^{do}$, with

[TABLE]

Here we have introduced the compact notation

[TABLE]

The diagonal matrix $\Delta_{d}$ contains the excitation energies $\Delta_{d,l}$ of $d_{r,\alpha,l}$ bosons. The $N\times N$ hopping matrix $M^{\Gamma/\Gamma^{\prime}}$ is a generalization of the fermion hopping parameter $T^{\Gamma,\Gamma^{\prime}}$, where $[M^{\Gamma/\Gamma^{\prime},\alpha\beta}_{-\bm{R}}]_{m,n}$ is the matrix element of $\Gamma/\Gamma^{\prime}$ interaction between states $\langle\Omega|d_{\bm{r},\alpha,m}$ and $d_{\bm{r}+\bm{R},\beta,n}^{\dagger}|\Omega\rangle$. E.g.,

[TABLE]

[TABLE]

The $1\times N$ tunneling vector $\mathcal{T}^{\alpha}$ is defined as the matrix element of $\Gamma^{\prime}$ interaction between NN and NNN boson states. E.g.,

[TABLE]

[TABLE]

Similar to the fermion vison pairs, the $\Gamma$-induced boson hopping is much stronger in the FM Kitaev model than in the AFM Kitaev model, making the boson band gap closes at much smaller $\Gamma$ values in the FM case. On the other hand, the $\Gamma^{\prime}$-induced boson hopping is comparable in both cases.

Once the boson gap closes at certain momentum $\bm{k}$, the resulting condensed boson mode may induce a specific magnetic long-range order, determined by $\bm{k}$ and the wave function of the soft mode—particularly its component on NN vison pairs. As the NNN vison pairs cannot be created or annihilated by local Pauli operators, their condensation does not immediately lead to appearance of the expectation value of local spin operators. However, because they can be generated by the $\Gamma$ interaction, their condensation affects the bond energies and may induce bond anisotropy (nematicity), particularly if their wavefunctions carry finite “orbital” angular momenta. We will show that this is likely the case when the FM KSL is destroyed by a positive $\Gamma$.

III FM Kitaev model

III.1 Consequences of the $\Gamma$ interaction

As shown in earlier sections, for a FM Kitaev interaction, the $\Gamma$ term induces strong single-vison hopping, causing the vison gap to close at relatively small $\Gamma$ values. Interestingly, we have found that the gap closes at an IC wave vector $\bm{k}\approx 0.67\,\mathrm{M}$ between the $\Gamma$ and $\mathrm{M}$ points when $\Gamma\approx 0.035$, and at the $\Gamma$ and $\mathrm{K}$ points when $\Gamma\approx-0.02$ (see Fig. 3 and Fig. 7(a)). These vison gap closings signal phase transitions into new phases driven by vison proliferation. However, the nature of the emergent phases is difficult to determine directly. This difficulty stems from several factors. First, since the model is gapless, the vison is a non-local excitation, making it unclear how to assign a self-statistical phase to it. Moreover, in nearby phases where the matter fermions are gapped by a Haldane mass, the vison becomes a local excitation but acquires a non-Abelian character. A microscopic understanding of what phase transitions can be driven by the softening of such non-Abelian particles remains a major open problem (see, however, Ref. [58] for a field-theoretic proposal of a direct continuous transition from non-Abelian states to short-range entangled states). We expect, however, that generally the closure of the single-vison gap drives confinement of all non-local quasiparticles, leading to the disappearance of topological order. This can be argued by noting that the vison “core” is a boson dressed by the non-local cloud of itinerant fermions [37, 59], and the condensation of this core is responsible for confinement, as it can be more clearly understood in Abelian cases [36]. This principle is supported by our numerical estimates for the critical $\Gamma$ values at which this occurs, which agree with previous numerical results for phase transitions into other ordered states [46, 21].

We have also found that the fermionic vison pairs acquire significant $\Gamma$-induced hopping, leading to a band gap closing at relatively small $\Gamma$ values. For positive $\Gamma$, the lowest $\tilde{\chi}$ band is flat and becomes gapless around $\Gamma\approx 0.05$. For negative $\Gamma$, the $\tilde{\chi}$ gap closes at $\Gamma$ and $\mathrm{K}$ points when $\Gamma\approx-0.09$ (see Figs. 5(a)–(b) and Fig. 7(a)). Notably, in the absence of $\Gamma^{\prime}$ interaction, the $\tilde{\chi}$ and matter $c$ fermions are decoupled, so the $c$-fermion Dirac cones persist and remain gapless along the $\Gamma$ axis.

The bosonic vison pairs also develop strong $\Gamma$-induced hopping, leading to band gap closures at critical $\Gamma$ values comparable to those of single visons and fermions. In the absence of $\Gamma^{\prime}$ interaction, the NN ($d$) and NNN ($o$) vison pairs are decoupled. Along the $\Gamma$ axis, the lowest-energy band originates purely from $o$ bosons (see Figs. 6(a)–(b)), whose gap closing does not immediately trigger local magnetic order, as these local bosons cannot be created or annihilated by single Pauli operators. This is also reflected in the absence of magnon-like low-energy dispersion in $S(\omega,\bm{Q})$ near the transition (see Figs. 9(a)–(b)).

For positive $\Gamma$, the lowest $o$ band is flat (becoming dispersive when $\Gamma^{\prime}$ is finite) and its gap closes around $\Gamma\approx 0.06$. There exists a twofold degeneracy at the $\Gamma$ point, with both modes carrying orbital angular momentum $l=\pm 1$:

[TABLE]

Because $o_{\bm{r},\alpha}$ boson can be created by the $\Gamma$ interaction on bond $(\bm{r},\alpha)$, i.e., $\langle\Omega|(\sigma_{\bm{r}}^{\beta}\sigma_{\bm{r}+\bm{\delta}_{\alpha}}^{\gamma}+\sigma_{\bm{r}}^{\gamma}\sigma_{\bm{r}+\bm{\delta}_{\alpha}}^{\beta})o_{\bm{r},\alpha}^{\dagger}|\Omega\rangle\neq 0$, the condensation of $\beta_{\Gamma,\pm 1}$ mode produces an additional contribution to the bond energy $(\Delta E^{x},\Delta E^{y},\Delta E^{z})\propto\Gamma(\cos(\theta),\cos(\theta+2\pi/3),\cos(\theta+4\pi/3))$. Depending on the Bose condensation phase $\theta$, $E^{\alpha}<E^{\beta}=E^{\gamma}$ or $E^{\alpha}>E^{\beta}=E^{\gamma}$ might occur. Interestingly, such nematic features were indeed observed in numerical studies of the $\Gamma$-induced paramagnetic (NP) phase [21, 47]. For $\Gamma<0$, the lowest $o$ band becomes dispersive and its gap closes at the $\Gamma$ and $\mathrm{K}$ points around $\Gamma\approx-0.04$. In this case, the soft mode at the $\Gamma$ point carries zero orbital angular momentum and does not induce bond-energy nematicity.

When $\Gamma^{\prime}$ is finite, the soft mode acquires weight on the NN $d$ bosons, and its condensation can induce long-range order with an average expectation value of single Pauli operators. The specific pattern of this order is governed by the momentum and wavefunction of the soft mode. In the phase diagram Fig. 7, the colored phase boundaries indicate the closure of the boson gap, with the color denoting the momentum of the soft boson mode: red for $\Gamma$, blue for $\mathrm{K}$, green for $\mathrm{M}$, and yellow for IC wave vectors. Near the $\Gamma$ axis at finite $\Gamma^{\prime}$, we find that the lowest-energy mode occurs at the $\mathrm{K}$ point, whose condensation yields a $120^{\circ}$ order (see its pattern in Fig. 8(d)). After applying the six-site $U_{6}$ transformation [50], this maps to a uniform FM order, and is therefore referred to as FM${}_{U_{6}}$ [14, 46, 32]. For negative $\Gamma$, it is known that the KSL should be replaced by an AFM${}_{U_{6}}$ order, which corresponds to the Néel state after the $U_{6}$ transformation [50, 14, 46]. While our analysis finds a lowest-energy boson mode leading to FM${}_{U_{6}}$ order, the second-lowest mode at $\mathrm{K}$ carries the appropriate structure to induce AFM${}_{U_{6}}$ order. We thus expect that improvements such as including bosonic states with two-$c$ fermion excitations for NNN vison pairs, or incorporating higher-order vison-pair hopping processes, might be needed to capture the correct magnetic phase. Note that without including the NNN vison pairs, the boson gap closes later than the fermion gap [40]. This highlights the important role of NNN vison pairs—as the lowest-energy two-vison configuration of the Kitaev model—in the low-energy physics.

III.2 Consequences of the $\Gamma^{\prime}$ interaction

To leading order, single visons exhibit vanishing $\Gamma^{\prime}$-induced hopping and thus remain gapped as $\Gamma^{\prime}$ increases, although higher-order terms may eventually enable their hopping and drive gap closure at certain critical $\Gamma^{\prime}$ values. While the $\tilde{\chi}$ and $c$ fermions become coupled in the presence of $\Gamma^{\prime}$ interaction, they remain decoupled at the $\mathrm{K}$ point, preserving the Dirac cone of the $c$ fermions. The fermion gap collapses at the $\Gamma$ point near $\Gamma^{\prime}\approx 0.24$, and at an IC wave vector around $\Gamma^{\prime}\approx-0.26$ (see Figs. 5(c)–(d)).

The boson band minimum lies at the $\Gamma$ point for both positive and negative $\Gamma^{\prime}$. Around $\Gamma^{\prime}\approx 0.17$, we find a FM order induced by boson condensation (see its pattern in Fig. 8(c) and $S(\omega,\bm{Q})$ in Fig. 9(d)), consistent with variational Monte Carlo results ($\Gamma^{\prime}_{c}\approx 0.12$ in Ref. [49]). For $\Gamma^{\prime}<0$, previous numerical studies have shown that the KSL is replaced by a zigzag order at $\Gamma^{\prime}\approx-0.1$ [20, 21, 47, 49], while our boson band calculations (see Fig. 6(c)) predict a gap closing near $\Gamma^{\prime}\approx-0.12$. Although the boson mode at the $\mathrm{M}$ point—corresponding to the zigzag wave vector—is not the lowest in energy (being higher than the $\Gamma$-mode that induces an out-of-plane FM order shown in Fig. 8(a)), its condensation leads to a magnetic pattern that precisely matches the zigzag order (see Fig. 8(b)). Furthermore, we find that the induced moment tilts at an angle of about $45^{\circ}$ from the honeycomb plane, close to the experimental value ($35^{\circ}$) observed in $\alpha$-RuCl3 [60]. These results suggest that bosonic vison pair condensation at the $\mathrm{M}$ point is the mechanism for the appearance of zigzag order in Kitaev materials, although further refinement—such as including $c$-fermion excited states for NNN vison pairs—is likely needed to to more accurately capture the correct energetics near the transition.

IV AFM Kitaev model

IV.1 Consequences of the $\Gamma$ interaction

In the AFM Kitaev model, single visons do not acquire $\Gamma$-induced hopping and thus remain gapped as $\Gamma$ increases. The fermions also have a weaker hopping, and their gap closes at much larger $\Gamma$ values (around $\Gamma\approx 0.38$ and $\Gamma\approx-0.8$) compared to the FM case.

On the other hand, the boson gap closes more rapidly than those of single visons and fermions. Unlike the FM case, where the lowest band along the $\Gamma$ axis originates from NNN vison pairs, the lowest band here is purely from NN vison pairs, whose condensation leads directly to finite expectation values of the local Pauli operators with magnetic long-range order. At $\Gamma\approx 0.28$, the boson gap closes simultaneously at the $\Gamma$ and $\mathrm{K}$ points (see Fig. 6(f)). We find that the soft mode at $\mathrm{K}$ point induces an in-plane FM${}_{U_{6}}$ order, as shown in Fig. 8(e), while the $\Gamma$ mode produces an out-of-plane Néel order (see Fig. 8(f)), which maps to an out-of-plane FM order under the $U_{6}$ transformation. Therefore, our results are fully consistent with the expected FM${}_{U_{6}}$ order [50, 46, 51, 32]. The dynamical structure factor $S(\omega,\bm{Q})$ near the transition point, shown in Fig. 9(f), exhibits a magnon-like low-energy dispersion touching zero at both the $\mathrm{K}$ and $\Gamma^{\prime}$ points. Notably, it closely resembles results from recent random-phase approximation calculations [61]. In the presence of $\Gamma^{\prime}$ interaction, the $\Gamma$–$\mathrm{K}$ degeneracy is lifted: a $\Gamma$-mode-driven Néel order emerges for $\Gamma^{\prime}>0$, while a $\mathrm{K}$-mode-driven FM${}_{U_{6}}$ order appears for $\Gamma^{\prime}<0$ (see Fig. 7(b)). Interestingly, a similar trend has also been reported in a vison-pair analysis of the $K$–$J$–$\Gamma$ model, with $\Gamma\textquoteright$ replaced by $J$ [40].

Around $\Gamma\approx-0.4$, we find that the boson band touches zero at an IC wave vector $\bm{Q}_{\mathrm{IC}}$ located between the $\Gamma$ and $\mathrm{M}$ points (see Fig. 6(e)), therefore the condensation of this mode is expected to induce an IC spiral order (see its pattern in Fig. 8(g) and dynamical structure factor in Fig. 9(e)). Moreover, a nearly gapless mode also appears at the $\mathrm{M}$ point, whose condensation leads to a stripy order, as shown in Fig. 8(h). The close proximity in energy between this mode and the IC mode indicates a strong competition between the two orders—consistent with recent numerical studies, where both types of order emerge depending on system size, and are collectively referred to as a “spin nematic” phase [46]. Upon introducing a small positive $\Gamma^{\prime}$, the $\mathrm{M}$ mode softens first, favoring a stripy phase; conversely, for negative $\Gamma^{\prime}$, the IC mode becomes soft first, leading to an IC spiral order where the value of $\bm{Q}_{\mathrm{IC}}$ evolves with $\Gamma^{\prime}$. In agreement with our findings, the IC order has also been obtained in a vison-pair analysis of the $K$-$J$-$\Gamma$ model with AFM Kitaev coupling and $\Gamma<0$ [40]. Morevoer, previous numerical studies have reported analogous behavior when $\Gamma^{\prime}$ is replaced by a Heisenberg coupling: the system exhibits stripy order for $J>0$ and IC spiral order $J<0$ [14].

It is worth noting that the condensation of bosonic vison pairs does not immediately lead to spinon confinement, as they are local bosons whose statistical interactions with matter fermions are trivial. Therefore, since we have found that single visons remain gapped due to the absence of $\Gamma$-induced hopping, the KSL is likely to evolve into a magnetically ordered QSL before undergoing a transition to a conventional ordered phase at larger $\Gamma$ values, where single-vison condensation occurs.

To further investigate this, we analyze the $\Gamma$-induced single-vison hopping at second order. We find that it also vanishes exactly in the AFM Kitaev model, suggesting that single visons are unlikely to condense via enhancing its mobility. However, we have found that at third order, a $\Gamma$-induced chemical potential term arises: $\propto-\Gamma^{3}\sum_{i}n_{i}^{\mathrm{v}}$ (see details in the Supplemental Material). This term favors (disfavors) single-vison excitations for $\Gamma>0$ ($\Gamma<0$). Including this term leads to single-vison gap closing for $\Gamma>0$, while they remain gapped for $\Gamma<0$ (see the gray dashed line in Fig. 7(b)). This analysis provides additional evidence that the FM${}_{U_{6}}$ order observed for $\Gamma>0$ corresponds to a conventional (confined) magnetic order, whereas the IC phase at $\Gamma<0$ may represent a magnetically ordered yet genuine QSL. Although no direct numerical evidence for the QSL has been reported, previous DMRG studies on a two-leg ladder identified this regime as a line of first-order transitions between stripy and IC phases, leaving open the possibility that such an exotic QSL could emerge in the two-dimensional limit [46].

IV.2 Consequences of the $\Gamma^{\prime}$ interaction

Different from the $\Gamma$ interaction, the $\Gamma^{\prime}$ interaction induces a finite single-vison hopping. Similar to the case of a Zeeman field, the $\Gamma^{\prime}$-induced single vison hopping has a fractionalized translational symmetry, resulting in two single-vison bands with opposite Chern number [38, 52]. The single-vison gap closes around $\Gamma^{\prime}\approx\pm 0.21$. The fermionic vison bands are also dispersive (see Fig. 5), with the gap closing first at the $\mathrm{M}$ point for $\Gamma^{\prime}\approx 0.1$ and at the $\Gamma$ point for $\Gamma^{\prime}\approx-0.1$.

At $\Gamma^{\prime}\approx 0.09$ and $\Gamma^{\prime}\approx-0.13$, the bosonic vison pairs close gap at the $\Gamma$ point. The condensation of the soft boson results in an out-of-plane Néel order, similar to the one shown in Fig. 8(f). The corresponding dynamical structure factors near both transitions exhibit similar features —a magnon-like low-energy band and a mushroom-shaped high-energy continuum (see Figs. 9(g) and (h))—which could be examined in future numerical studies.

V Discussion

We have investigated the $K$-$\Gamma$-$\Gamma^{\prime}$ model from the perspective of the quasiparticles of the KSL, treating the non-Kitaev $\Gamma$ and $\Gamma^{\prime}$ interactions as perturbations and analyzing their effects on the quasiparticles’ dispersions. This quasiparticle-based approach has previously been applied to other extended Kitaev models, including the $K$-$J$-$\Gamma$ model [40], Kitaev model with external magnetic fields [52, 38, 41], and multilayer Kitaev systems [62], offering valuable insights into their physical properties and the competitions between their QSLs and other magnetically ordered phases.

The stronger (weaker) $\Gamma$-induced hopping of single visons and vison pairs in the FM (AFM) Kitaev model naturally explains the weaker (stronger) stability of the KSL under $\Gamma$ perturbations, as reflected in the phase diagram shown in Fig. 7. This suggests that the suppression of the KSL is driven by vison proliferation.

The softening of bosonic vison pairs provides a compelling mechanism for the emergence of various magnetic orders, which we have been able to infer by analyzing the momentum and other quantum numbers of the wavefunction of softened bosonic modes. For example, in the AFM Kitaev model, the FM${}_{U_{6}}$ order at $\Gamma>0$ arises from the condensation of the $\mathrm{K}$ mode or $\Gamma$ mode (see Figs. 8(e)–(f) and Fig. 9(f)); the IC spiral order and the stripy orders at $\Gamma<0$ are driven by softened modes at $\bm{Q}_{\mathrm{IC}}$ and $\mathrm{M}$ points (see Figs. 8(g)–(h) and Fig. 9(e)), respectively. In the FM Kitaev model with $\Gamma^{\prime}<0$, although the $\mathrm{M}$ mode—which induces a zigzag order consistent with the moment direction observed in $\alpha$-RuCl3—is only the second lowest-energy mode (see Fig. 6(c) and Fig. 8(b)), its energy might be lowered by extending the vison-pair subspace, possibly allowing it to become the softest mode.

The NNN vison pairs ($o$) also play an essential role at low energies. In the FM Kitaev model, they form the lowest bosonic band along the $\Gamma$ axis in the $\Gamma$-$\Gamma^{\prime}$ phase diagram (see Fig. 7(a)). For $\Gamma>0$, the softened mode at $\Gamma$ point carries finite orbital angular momentum and induces a paramagnetic phase with anisotropic bond energies, offering an explanation for the nematic paramagnetic phase observed in previous numerical studies [21, 47].

One of our main finding is that, in the AFM Kitaev model with $\Gamma<0$, single visons remain gapped and matter fermions stay deconfined even after the bosonic vison pairs condense (see Fig. 7(b)). This raises the tantalizing possibility of an intermediate magnetically ordered QSL before the system transitions into a trivial magnetically ordered phase at larger $\Gamma$ values. Our findings therefore allow us to conjecture that the AFM Kitaev model with $\Gamma\lesssim-0.4|K|$ realizes a QSL that retains the anyon content of the KSL while simultaneously exhibiting broken symmetry. The symmetry-breaking magnetic order is likely either stripy or a IC spiral (see Figs. 8(g)–(h)). Consistently, the IC order has also been reported in a previous vison-pair study of the $K$-$J$-$\Gamma$ model [40]. The coexistence of fractionalization and a spiraling magnetic moment is particularly intriguing, as the IC moment may remain strictly gapless. We hope that future numerical and analytical studies will provide a more direct assessment of this fascinating possibility.

We have also demonstrated how correlation functions of local observables can be directly related to the bosonic vison pairs, and used these modes to explicitly compute the dynamical structure factor. Near phase transitions into magnetically ordered states, the dynamical structure factor exhibits low-energy magnon-like dispersions along with broad high-energy continua. Notably, in some cases, our results resemble those obtained from random-phase approximation studies (cf. Fig. 9(f) and Fig. 9(c) of Ref. [61]). Such features serve as valuable signatures for probing emergent phases in Kitaev materials. This framework can also be extended to study other experimental probes, such as Raman spectroscopy [63], providing a promising direction for future investigations.

VI Methods

VI.1 Magnetic orders from the condensation of bosonic vison pairs

The condensation of a bosonic mode can give rise to magnetic long-range order, provided it carries weight on NN vison pairs (see also Sec. V of Ref. [40]). Consider a soft boson mode $\beta_{\bm{q}}$ with momentum $\bm{q}$:

[TABLE]

where the $N\times 1$ vector $w_{\alpha}(\bm{q})$ is its weight on NN $d_{\bm{q},\alpha}$ vison pairs and $v_{\alpha}(q)$ is its weight on NNN $o_{\bm{q},\alpha}$ vison pairs. Note $d_{\bm{q},\alpha}^{\dagger}=(d_{\bm{q},\alpha,1}^{\dagger},\dots,d_{\bm{q},\alpha,N}^{\dagger})$ (see Eq. (13)), and

[TABLE]

The local magnetic moment induced by the condensation of $\beta_{\bm{q}}$ mode reads:

[TABLE]

Here $\theta$ is the condensation phase of the soft mode. It can be shown that, for $\bm{r}\in A$,

[TABLE]

Here

[TABLE]

Similarly, for $\bm{r}^{\prime}\in B$, there is

[TABLE]

with

[TABLE]

Therefore, once $w_{\alpha}(\bm{q})$ and $P_{A/B,\alpha}$ are determined, the resulting local moments can be obtained. In our calculations, we find that for orders induced by the $\Gamma$ or $\mathrm{M}$ modes—such as FM, Néel, zigzag, and stripy—the direction of the local moments is independent of $\theta$. In contrast, for the $120^{\circ}$ order induced by the $\mathrm{K}$ mode, the local moment direction depends on $\theta$, consistent with Ref. [40].

VI.2 Dynamical structure factor from the bosonic vison pairs

In this section, we discuss the dynamical structure factor approximated from the bosonic vison pairs, a method first used in Ref. [40]. The dynamical structure factor is defined as

[TABLE]

where $N_{s}=2N$ denotes the total number of sites in the system. The Fourier tranform of Pauli operators is given by

[TABLE]

Note that due to the bipartitle nature of the honeycomb lattice, $\sigma_{\bm{Q}}^{\alpha}\neq\sigma_{\bm{Q}+\bm{b}}^{\alpha}$ with $\bm{b}$ being a reciprocal lattice vector. As a result, it is helpful to evaluate $S(\omega,\bm{Q})$ in the extended BZ.

$|\mathrm{GS}\rangle=|\Omega\rangle$ before the breakdown of KSL. Using the Lehmann representation, Eq. (VI.2) can be recast into:

[TABLE]

with $E_{l}-E_{0}$ being the excitation energy of $l$-th excited state. Since a local Pauli operator always excites a single bosonic vison pair from the ground state $|\Omega\rangle$ [54, 55], one can choose $|l\rangle$ to be the single-boson eigenstate $\beta_{\bm{Q},l}^{\dagger}|\Omega\rangle$, where $\beta_{\bm{Q},l}^{\dagger}$ creates the $l$-th eigenmode with momentum $\bm{Q}$ of the boson Hamiltonian. Note that for $\bm{Q}$ outside the first BZ, one should use $\beta_{[\bm{Q}],l}^{\dagger}|\Omega\rangle$, where $[\bm{Q}]$ denotes the corresponding wave vector folded back into the first BZ. After some algebra, one can show that:

[TABLE]

The $S(\omega,\bm{Q})$ thus reads:

[TABLE]

The Dirac delta function is approximated by a Lorentzian: $\delta(x)\approx\frac{\eta/\pi}{\omega^{2}+\eta^{2}}$, where we take $\eta=0.04|K|$ in our calculation. Fig. 10 presents the dynamical structure factor $S(\omega,\bm{Q})$ of the ideal Kitaev model, computed from bosonic vison pairs. The results show good agreement with the exact solution of the Kitaev model [54, 55].

Acknowledgements.

We thank Chao-Ming Jian, Aprem Joy, Peng Rao, and Jiucai Wang for helpful discussions. C.C. is grateful to Shang-Shun Zhang for several insightful discussions. C.C. acknowledges support from the National Natural Science Foundation of China (Grants No. 12404175 and No. 12247101), the Fundamental Research Funds for the Central Universities (Grant No. lzujbky-2024-jdzx06), the Natural Science Foundation of Gansu Province (No. 22JR5RA389), and the ‘111 Center’ under Grant No. B20063. I.S.V. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) under research Grants No. 542614019, No. 518372354, and No. 555335098. We thank Beijing Paratera Co., Ltd. for providing HPC resources that contributed to the numerical results reported in this paper.

Data availability

The datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.

Code availability

The codes used during the current study will be available from the corresponding author upon reasonable request.

Author contributions

C.C. and I.S.V. designed the study. C.C. carried out the calculations with input from I.S.V. C.C. and I.S.V. wrote the manuscript together.