Kinetic Ferromagnetism and Topological Magnons of the Hole-Doped Kitaev Spin Liquid

We study the effect of hole doping on the Kitaev spin liquid (KSL) and find that for ferromagnetic (FM) Kitaev exchange $K$ the system is very susceptible to the formation of a FM spin polarization. Through density matrix renormalization group (DMRG) simulations on finite systems, we uncover that the introduction of a single hole with a hopping strength of just $t\sim{}0.28K$ is enough to disrupt fractionalization and polarize the spins in the [001] direction due to an order-by-disorder mechanism. Taking into account a material relevant FM anisotropic spin exchange $\varGamma$ drives the polarization towards the [111] direction via a reorientation transition into a topological FM state with chiral magnon excitations. We develop a parton mean-field theory incorporating fermionic holons and bosonic spinons/magnons, which accounts for the doping induced FM phases and topological magnon excitations. We discuss experimental signatures and implications for Kitaev candidate materials.


INTRODUCTION
Kinetic ferromagnetism, resulting from the subtle interplay between the motion of electrons and their interactions, provides a counter-intuitive example of a strong interaction effect in condensed matter physics [1,2].Nagaoka [3] famously showed how the interference of paths from a single hole doped into a half-filled Hubbard model with infinite repulsion U can lead to a ferromagnetic (FM) state in a system that typically supports antiferromagnetic (AFM) order.Given that the required large interaction limit is an experimental challenge, signatures of Nagaosa's ferromagnetism have only recently been observed experimentally in quantum dots [4] and semiconductor heterostructures [5].In this context, intriguing questions are whether similar kinetic magnetism can appear in other correlated models that are experimentally accessible; and whether the kinetic FM state itself can be nontrivial, e.g., host chiral excitations.
In this work, we account for the effects of hole doping by considering an extension of the Kitaev honeycomb model, referred to as the t-K model, where t denotes the hopping strength of doped holes.In the slow hole limit (t ≪ K), the KSL phase is robust and holes only introduce quasi-static vacancies [49,50], see Fig. 1(a).However, in realistic systems, t is typically much larger than the Kitaev exchange K, and previous works based on parton meanfield theories suggested superconducting ground states in the hole-doped regime [51][52][53][54].Yet, our numerics does not support superconductivity at the considered low doping values, see Supplementary Note 2. Recently, a density matrix renormalization group (DMRG) study suggested a charge density wave ground state emerging on top of a hole-doped AFM KSL, in which the superconducting correlations fall off almost exponentially at long distances [55].
We investigate the ground state of the t-K model by DMRG [56,57] and show that the FM KSL is remarkably fragile already for small and slow hole doping.For hopping strengths t of the order of the KSL's flux gap, ∆E v ∼ 0.1K, the system is already partially FM polarized by the itinerant holes, spontaneously breaking the time-reversal symmetry, as illustrated in Fig. 1(b).To account for our numerical results, we develop a parton theory incorporating fermionic holes and bosonic spinons/magnons.It unveils that the hole kinetic term effectively serves as a FM Heisenberg coupling destabilizing the KSL.In addition, our parton theory shows that the resulting FM order along the [001] direction originates from an order-by-disorder mechanism [58,59].We further uncover that the presence of a finite FM off-diagonal exchange, Γ > 0, shifts the magnetization direction from [001] to [111].Remarkably, due to the change of spin polarization direction, the kinetic FM state becomes non-trivial.It spontaneously forms topological FM order with chiral magnon edge modes, akin to the FM state of the Kitaev model induced by a strong external magnetic field [60,61].Furthermore, we show that hole doping can significantly lower the critical field above which the field-polarized FM state appears in Kitaev candidate materials, such as α-RuCl 3 .

Model Hamiltonian
The celebrated Kitaev honeycomb model is defined by the following Hamiltonian where S γ j (γ = x, y, z) are the S = 1/2 spin operators and S γ j S γ k are Ising couplings according to the γ-type of nearestneighbor (NN) ⟨ jk⟩ bonds.Microscopic derivations [32,33,[62][63][64] have shown that the Kitaev interaction of several material candidates is likely FM.Therefore, we focus on positive K and set K = 1 as the unit of energy.There exist conserved plaquette operators W p on each hexagon p as where the lattice sites j = 1, ..., 6 correspond to the elementary hexagon plaquette p shown in Fig. 1(a).The pristine KSL ground state has W p = 1 for all plaquettes.
To describe the physics of hole doping, we introduce the t-K model [51,65,66] where c † s is the creation operator of an electron with spin index s =↑, ↓ and the projector P removes doubly occupied states.
Spin operators, S γ j = c † j σ γ c j /2, are given by the fermionic vectors c † j = (c † j,↑ , c † j,↓ ) and the standard σ γ Pauli matrices.The hole doping is parameterized by δ such that s ⟨c † j,s c j,s ⟩ = 1−δ.Note that the NN hopping is spin-independent and the effects of spin-orbital coupling are nevertheless retained in H S [65], resulting in a space symmetry group of D 3d .Besides, the t-K model is also symmetric under time reversal transformation and preserves the charge-U(1) conservation.
In order to implement DMRG calculations, the system is placed on a two-dimensional cylindrical geometry, dubbed YCL y , with periodic boundary conditions (PBCs) along the short direction (with L y unit cells), while the longer one (L x ) is open.See details about the honeycomb lattice in Supplemental Note 1.

Phase Diagram
First, we focus on the one-hole doped system.The groundstate phase diagram on YC4 cylinders is analyzed by studying the average magnetic moments = ⟨S i ⟩ and the average Z 2 flux W p , as shown in Fig. 1(c,d).For t < t * , the ground state corresponds to a site-diluted KSL, where the kinetic energy of the holes is insufficient to overcome the Z 2 flux excitation gap ∆E v ∼ 0.1K.Instead, the presence of the slow-moving holes can be thought of as quasi-static vacancies within the KSL state [49,50].Although W p slightly decreases as t increases in this phase, the system does not acquire any magnetization.
We observe that a phase transition occurs at t = t * ≈ 0.28K, as visible by kinks in the curves of both the magnetization and Z 2 flux, see inset plots of Figs.1(c) and (d).As one of our main findings, already for t > t * a FM phase appears with a finite magnetization along the [001] direction.We note that the first-order derivative of the ground-state energy also exhibits a kink around t * ≈ 0.28K which suggests that the phase transition is of second order, see Supplementary Note 2.
This FM phase is an example of kinetic magnetism, which spontaneously breaks the time-reversal symmetry and arises solely from the holes' kinetic energy.Remarkably, the KSL correlations are only partially depleted, as signaled by the gradual decrease of the Z 2 flux W p in Fig. 1(d).This coexistence of magnetization and finite plaquette flux is a general feature of the [001] FM phase until the system is fully polarized by a sufficiently large hopping strength.According to Nagaoka's theorem, for a single-hole doped system [3] in the thermodynamic limit, the saturated FM state is only reached for vanishing spin exchange (or in other words, the on-site Hubbard U → ∞).Therefore, the intriguing coexistence of FM order and KSL correlations could generally persist in systems described by the t-K model.
We further investigate whether FM order persists for multiple holes and study the magnetization as a function of hole density δ in Fig. 2. Unlike Nagaoka's ferromagnetism, which may disappear for a thermodynamic density of holes [67], we find that the FM order observed in our work is robust and can be further enhanced with multi-hole doping.First, a larger δ can lower the critical value of t * .For instance, at t = 0.2K, the system remains in the site-diluted KSL phase without any magnetization for δ ≲ 0.02, while a small but finite magnetic order emerges for δ > 0.02.For the hopping strength t which is already capable of polarizing spins at δ ≈ 0.01, a slight increase in δ can significantly enhance the magnetization until it reaches a saturation value at δ ≈ 0.05.Due to a rapid increase in entanglement entropy, our DMRG simulations face convergence problems for systems with even higher doping levels.Nevertheless, our results suggest that the ferromagnetic order is a generic feature of the t-K model, at least in the low doping limit δ ≤ 0.06.

Off-Diagonal Exchanges
Next, we study the hole-doped extended Kitaev honeycomb model, dubbed the t-K-Γ model, as a more realistic system that incorporates additional off-diagonal symmetric spin exchange parameterized by Γ.The spin part of the model now becomes with α β γ for the Γ terms on γ-type bonds.Although the Γ term is typically believed to be AFM [46,68], here we focus on the effect of a FM Γ > 0 as this has the most interesting consequences.
In the absence of Γ, the spin polarization is always along the [001] direction (or by symmetry equivalently, the [010] and [100] directions), which can be attributed to an order-bydisorder mechanism as we will explain below.In Fig. 3(a), we show that the presence of Γ > 0 can shift the spin polarization direction from [001] to [111].Indeed, already a tiny value of Γ * is sufficient to induce a significant change in the polar angle of magnetization, denoted as ϕ ≡ cos −1 (M z /|M|).For Γ = 0.02, the polar angle is already very close to the value of , which corresponds to a FM order along the [111] direction, see Fig. 3(b).At Γ * ≈ 0.01K, a noticeable kink in the first-order derivative of the ground-state energy in Fig. 3(c) indicates that this change in spin polarization direction is accompanied by a phase transition.We will discuss below that this transition leads to topologically nontrivial magnon excitations.The magnetization displays a dip around Γ * , corroborating the findings from the analysis of the ground-state energy.Note that a moderately large Γ can also enhance the magnitude of the magnetization.
Compared to the FM Γ, the perturbation of an AFM Γ has a much smaller impact on the [001] ordered ground state.For instance, a small AFM Γ ≈ − 0.02 only induces a canting field in the [001] order.However, a FM Γ with a similar magnitude is sufficient to shift the polarization direction from [001] to [111].

Parton mean-field theory
In order to understand the selection of the FM moment direction and the impact of Γ, we develop a parton meanfield theory.To effectively describe the low-energy degrees of freedom in the kinetic FM phase, we consider the following parton representation for electron operators [69] with c † j,s ≡ h j b † j,s .Here, b † j,s 's are Schwinger boson operators representing spinons, and h j 's are fermionic holon operators representing empty sites.The local constraint, h † j h j + s b † j,s b j,s = 1, ensures the original three-dimensional physical local Hilbert space.The parton theory is constructed such that it allows for spontaneous symmetry breaking toward ordered states.
The t-K-Γ model takes the form with F jk ≡ s b † j,s b k,s representing short-range FM spin correlations [70].One can observe that the hopping term in Eq. ( 4) connects the hole's kinetic energy P jk ≡ h † j h k with FM correlations F jk , suggesting that finite doping with dominant t dramatically renormalizes the spin-spin interactions and fosters FM order.
The second term in Eq. ( 4), Ĥγ jk , represents spin interactions as defined in Eq. ( 3), but expressed in terms of bosonic spinons instead of electron operators.We develop a large-N mean field theory within a Schwinger-boson approach [70][71][72] in which the spinon occupancy is given by n b = s b † j,s b j,s = 1 − δ (see Supplementary Note 3).Guided by our DMRG results, our interest lies in the FM phase with t ≫ K.We also explored the possibility of a gapped QSL phase but found no evidence thereof in the limit of n b → 1.Hence, the most natural ansatz is a FM state along the general direction of (sin ϕ cos θ, sin ϕ sin θ, cos ϕ).For conveniently describing magnon excitations, we can replace the Schwinger bosons with Holstein-Primakoff (HP) bosons as [1,73] where a † j and a j are the creation and annihilation operators for the HP bosons (magnons), respectively.Then we perform a standard Hartree-Fock decoupling of Eq. ( 4) to obtain a meanfield theory with separate hole and spin parts, dubbed H h and H s , respectively.Here, H h refers to a spinless free-fermion band with renormalized bandwidth by the factor |⟨ F jk ⟩| and filling of δ.The spin part H s is treated in spin-wave theory (SWT) incorporating up to fourth-order Holstein-Primakoff expansions, and the spin amplitude is renormalized to S = n b /2.The expectation values of FM spin correlations F jk and hole's kinetic energy P jk are determined self-consistently.Remarkably, we can now see that for δ ≪ 1 the kinetic energy effectively acts as a FM Heisenberg interaction with a coupling constant J ≡ − ⟨P jk ⟩t/n b ∼ |t|δ.
The semi-classical ground-state energy of the FM Kitaev honeycomb model (t = Γ = 0) turns out to be independent of ϕ and θ, yielding an emergent O(3) manifold of degenerate ground states.This classical O(3) degeneracy can be lifted by quantum fluctuations in a quantum order-by-disorder mechanism [58,59].Within SWT the quantum correction to the ground-state energy is , where N c is the number of the unit cells and ω k,n are the two magnon bands of the honeycomb lattice.Note that ω k,n depend on ϕ and θ implicitly, and thereby so does E s .Our calculations reveal that a magnetization along the [001] direction (or equivalently, the [100] and [010] directions) is energetically preferred in accordance with our DMRG results (see Supplementary Note 3).
A complication arises from the fact that the semi-classical ground state of the pure Kitaev honeycomb model is a classical spin liquid [74] which is manifest in the linear SWT for arbitrary spin polarization directions as a nearly flat magnon band with almost zero energy (see Supplementary Note 3).Thus, an HP expansion that only includes the quadratic terms is insufficient to correctly capture the quantum fluctuations.We employ a nonlinear SWT by including the quartic terms in the HP expansion and treat these within the Hartree-Fock decomposition [75].We indeed find significant renormalizations of the flat magnon bands, for instance, a closing of the gap between formerly flat magnon bands in the nonlinear SWT for Γ = 0.More technical details can be found in Supplementary Note 3.
Including a FM Γ > 0 breaks the semi-classical O(3) degeneracy of the ground state.Our SWT shows that the semiclassical ground state has [111] magnetic order rather than [001], which is again consistent with our DMRG calculations, see Fig. 3(b).From the magnetic field polarized regime of the Kitaev model, it is known that different polarization directions can have qualitatively different types of boundary magnon excitations as first pointed out in Ref. [60,61].For spins polarized along the [001] direction with Γ = 0, the two magnon bands touch linearly and are almost non-dispersive along one of two primitive vectors.The system thereby cannot support nontrivial boundary excitations.On the other hand, the introduction of a finite Γ > 0 modifies the magnon bands from the change of the spin polarization to the [111] direction.Consequently, the two magnon bands are separated by a gap and acquire non-zero Chern numbers (−1, 1) for the lower and upper magnon bands, respectively [60,76,77].Hence, as our second main result, we find that the kinetic FM state of the doped Kitaev model (with small Γ > 0) spontaneously breaks TRS forming a topological magnon insulator.This explains the phase transition found when turning on the off-diagonal Γ interactions, see Fig. 3.
To see how the change in spin polarization can give rise to chiral boundary states, we calculate the magnon bands within SWT on narrow cylinders.With PBCs along y-direction, k y is still a good quantum number and the 2D systems can be treated as L y independent 1D chains.The magnon spectra for the [001] and [111] ordered states are shown in Figs.4(a) and (d).Indeed, only the [111] ordered state has chiral boundary modes (green) as predicted.

Real-time dynamics
So far our observation that a topological FM state is realized for the [111] direction solely relies on the parton+SWT description.
A numerical confirmation by DMRG is complicated by the fact that the ground state is topologically trivial but only the magnon excitation spectrum carries signatures of the magnon Chern bands.Therefore, we study the distinct real-time dynamics of excitations in the [001] and [111] ordered FM states.To obtain an initial state with a magnon-like excitation at the boundary, we apply a local time-reversal operator at lattice site j, which couples to spin excitations for arbitrary spin polarizations, to a one-hole doped ground state |Ψ⟩ as where K denotes the complex conjugate operator.Note that σ y j also annihilates the components of |Ψ⟩ which contain a local hole at site j.As a result, in addition to flipping the spins at site j, it also weakly excites charge excitations and slightly modifies the global magnetization pattern.Starting with the prepared initial state |Ψ( j)⟩, the system is evolved by performing a standard real-time evolution as |Ψ( j, t)⟩ = e −iHt |Ψ( j)⟩.In practice, we approximately represent the shorttime unitary operator, e −iHdt , as a compact matrix product operator (MPO) [78], where a small time step dt = 0.01K has been used.Then the time evolution can be efficiently simulated by applying such an MPO to |Ψ( j, t)⟩ successively.
We create a one magnon defect at site j in the leftmost column of the YC cylinders and follow the spin current for each column x defined as with ⟨..⟩ the expectation value with respect to |Ψ( j, t)⟩ and (x, y) denoting the lattice coordinates.Note that j s (x, t) is obtained by summing over all the sites belonging to the xth column, namely, summing over allowed k y .Additionally, we also keep track of the charge current where n h j = 1 − P s c † j,s c j,s P is the occupancy of holes.The simulations for the dynamics are performed on YC3 cylinders, and the results are shown in Fig. 4. We observe a distinct difference in the spin currents between the [001] and [111] ordered states, as evidenced by the leftmost boundary columns in Figs.4(b) and 4(e), respectively.The spin current in the [001] ordered state shows very little dynamics, as indicated by its small magnitude of j s ∼ 10 −4 in Fig. 4(b).It can be attributed to the absence of distinct boundary modes and the almost vanishing velocity of magnon modes.In stark contrast, for the [111] direction j s (x, t) exhibits a strong and long-live signal in the leftmost boundary as expected for a localized surface excitation.Note that Eq. ( 5) unavoidably involves excitations related to bulk magnon bands, which explains why the boundary spin current can leak into the bulk.We have checked that the magnetic field-induced topological FM in the [111] direction [60,61] shows a very similar response confirming the presence of chiral magnons.
The charge current for the [001] and [111] ordered states are very similar, see Figs. 4(c) and (f), respectively.This is consistent with our parton mean-field theory.Both the [001] and [111] ordered states give rise to similar short-range FM spin correlations (see Supplementary Note 3), which suggests that the holon Hamiltonians are similar.
Implications for the field-polarized state of α-RuCl 3 Unlike the simplified models introduced in Eq. ( 1) and Eq. ( 3), more realistic models for experimentally relevant materials are much more complicated.Here, we focus on α-RuCl 3 and related materials [31,33,79] which are often described by a number of parameters, including {K, J, Γ, Γ ′ , J 3 }.The corresponding Hamiltonian is expressed as where α β γ and ⟨i j⟩ 3 denotes the third NN bonds.10); see also model parameters in Eq. ( 9).The left inset shows a typical zigzag order at (t, h in ) = (0.3, 0.03) and the right one shows a polarized state at (t, h in ) = (0.5, 0.06).
Positive and negative local magnetizations, ⟨S x j ⟩, are represented by red and blue dots, respectively, with their radii scaled according to the magnitude of ⟨S x j ⟩, with ⟨S x j ⟩ ≈ 0.21.The doping level is δ = 0.0625.
we focus on a specific point in the parameter space, given by which is close to a set of parameters identified in Ref. [79] for α-RuCl 3 , resulting in a zigzag-ordered ground state obtained by DMRG simulations.Moreover, an in-plane magnetic field suppresses the long-ranged zigzag order and ultimately polarizes the state for sufficiently strong fields.Furthermore, it has been proposed that an intermediate KSL state may separate these phases.To account for this effect, we introduce an additional Zeeman term in Eq. ( 8).The total spin Hamiltonian is then given by where h in represents the in-plane magnetic field along the [1 10] direction.
We then investigate how hole doping modifies the phase boundary of the zigzag phase.We add the same kinetic energy term for the holes as in Eq. ( 2), to the experimentally relevant Hamiltonian Eq. ( 10) and compute the ground-state phase diagram with DMRG; see Fig. 5. Remarkably, we find that the phase boundary between the zigzag and polarized phases is highly sensitive to hole doping.For a vanishing hopping strength, t → 0, the critical magnetic field required to polarize the spins is approximately h * in ≈ 0.085.However, already for a small hopping t = 0.8K and low hole doping of δ = 0.0625, the critical magnetic field is reduced to h * in ≈ 0.04.A larger value of t ≈ 1.8K can drive the zigzag order to a FM order (not shown).This can be understood by noting that |t|δ effectively acts as a FM Heisenberg exchange constant J.The polarized phase is distinct from the kinetic FM phase, as the latter develops magnetic order spontaneously.The significant effect induced by the kinetic energy of holes, therefore, is a complicating factor in identifying to correct Hamiltonian parameters in comparison to experimental measurements.

DISCUSSION
In summary, we have studied the hole-doped Kitaev-Γ model using DMRG and effective parton descriptions.It is by now well known that the FM KSL is much more fragile with respect to the application of an external magnetic field compared to the AFM KSL [80].Similarly, our DMRG simulations suggest that the AFM KSL is also much more robust to hole doping.Meanwhile, distinct hole spectral functions have been numerically observed in FM and AFM KSLs, in which a dynamical Nagaoka ferromagnetism emerges in the FM KSL only [81].Concentrating on FM Kitaev exchange, we have shown here that small hole doping is sufficient to destabilize the KSL leading to a kinetic FM state generally coexisting with finite plaquette fluxes.We emphasize that, unlike the field-polarized phase in a magnetic field [82], the doping-induced FM order breaks the timereversal symmetry spontaneously.Due to the particle-hole symmetry of our model, our conclusions also hold for charge doping.
We developed a parton mean-field theory, incorporating fermionic holons and bosonic magnons, which indeed shows that hole condensation gives rise to effective FM Heisenberg exchanges.We found that in the absence of off-diagonal exchange Γ, the spin polarization is along the [001] direction due to an order-by-disorder mechanism, which is also supported by DMRG simulations.A FM Γ term switches the polarization direction to [111], which leads to a distinct topological FM phase with chiral boundary magnon excitations.
Our work raises a whole range of questions for future research.First, one could investigate flux threading and topological entanglement of the FM phase coexisting with finite plaquette fluxes to probe the presence of topological order.Second, it would be interesting to search for other kinetic FM phases of Hubbard models hosting chiral magnon excitations.Third, to further investigate the transition between the KSL and kinetic FM phases, it would be desirable to develop a parton theory incorporating bosonic holons and fermionic spinons.In this context, such a theory could be the starting point for a self-consistent random phase approximation [83] or Gutzwiller-boosted DMRG [84][85][86][87] in order to study the ordering instabilities and spin excitation spectra.Forth, in our real-time evolution results we could observe that the excitations of spin induce a (weak) response of charge, and vice versa.This could possibly be captured by a time-dependent parton mean-field theory in future studies to explore the interplay of charge and spin transport [88].Fifth, a study of doping effects for realistic model Hamiltonians could unveil a new mechanism (see Supplementary Note 4) for stabilizing the zigzag magnetic orders observed experimentally.This could also be important for refining the microscopic models of RuCl 3 and its peculiar magnetic field response, both of which remain poorly understood [79].Sixth, the origin of the observed thermal Hall effect of RuCl 3 [89] is hotly debated [90] and it would be very interesting to explore how the topological magnons of the weakly doped Kitaev model studied here may contribute.
In conclusion, the Kitaev honeycomb model continues to be a fertile soil for novel physics.We expect the interplay of Kitaev spin exchange with kinetic electron motion to hold further surprises in the future.

Supplemental Materials for "Kinetic Ferromagnetism and Topological Magnons of the Hole-Doped Kitaev Spin Liquid"
This Supplemental Material includes more details of the DMRG calculations and parton theory.In Supplementary Note 1, we show the considered honeycomb lattice and the corresponding one-dimensional path for DMRG simulations.In Supplementary Note 2, we show additional information about the DMRG results.In Supplementary Note 3, we show details of the parton mean-field theory.In Supplementary Note 4, we provide additional results for the zig-zag order.In Supplementary Note 5, we provide more information about the spin-1/2 moment around the exactly solvable point of t = 0. To employ the MPS-based methods, one must establish a site ordering for the honeycomb lattice.This can be achieved by assigning an integer from 1 to 2N c to each lattice site, where N c is the number of unit cells.As shown in Fig. 6, we utilize a site-labeling scheme for the honeycomb lattice on a YC4 (L y = 4) cylinder with length L x with primitive vectors ⃗ x and ⃗ y.Note that the x-, y-, and z-type bonds for a Kitaev honeycomb model are marked by the blue, green, and red lines, respectively.

Ground-state calculations
For the DMRG simulations we make explicit use of a U(1) quantum number to control the level of hole doping.The bond dimension of DMRG is kept as large as χ = 4000, resulting in a typical truncation error of ≈ 7×10 −7 .To analyze the convergence of our DMRG simulations, we study the scaling behavior of the ground-state energy and magnetization with the inverse DMRG bond dimension 1/χ.As shown in Fig. 7, the ground-state energy varies weakly with χ for χ > 1000, suggesting a faithful convergence in our DMRG simulations.In the slow-hole regime (Kitaev spin liquid phase), the magnetization fluctuates around the magnitude of truncation error, indicating the absence of magnetic order.In the fast-hole regime (FM-ordered phase), the magnetization extrapolates to a finite value when 1/χ → 0, indicating that the kinetic FM order is stable when increasing the bond dimension.
Fig. 8 displays the ground-state energy E and its first-order derivative ∂ t E as functions of hopping strength t.Note that ∂ t E in Fig. 8(b) exhibits a kink at t * ≈ 0.28K which suggests that the phase transition is of the second order.

A. Superconducting correlation functions
The multiple-hole doping allows us to analyze the equal-time pairing-pairing correlation function.A diagnostic of the superconducting order is the pair-pair correlation function, defined as where ∆ † a (r 0 ) is SC pair-field creation operator with pairing index a defined on the γ-type (γ = x, y, z) NN bond r 0 and r is the distance along the cylinder (the x direction).We consider one even-parity pairing as where i 0 and j 0 denote two lattice sites on the NN bond r 0 .Similarly, we define three odd-parity pairing functions as In practice, we set i 0 and j 0 being the reference bond at x ∼ L x /4 to minimize the boundary effect.We study the abovementioned pair-pair correlation functions on the x-, y-, and z-type bonds for δ ≈ 0.04 and as shown in Fig. 9, we find that all of  11)-( 13)] obtained on a YC4 cylinder with Γ = 0, δ ≈ 0.04, and χ = 4000.For (a-c) we set t = 20K and for (d-e) t = 2K.The even-parity correlation functions, which are not shown here, are significantly smaller than the odd-parity functions (Φ e ∼ 10 −12 ).them exhibit exponential decays, The correlation lengths χ a for different pairing types are summarized in Fig. 9.This observation suggests the absence of superconductivity in the weak doping limit.

Real-time dynamics
The real-time dynamics is performed by the standard MPO-evolution algorithm [78].A small time step dt = 0.01[K] −1 has been used and At each intermediate step of the real-time evolution, one needs to truncate the MPO-evolved MPS.In order to estimate the accuracy of the final MPS, we introduce the accumulated truncation error defined by where ϵ j (D) is the sum of the discarded squared singular values at the j-th bond.Meanwhile, we also track the growth of entanglement entropy (EE) of the evolved MPS.As shown in Fig. 10, the truncation errors are close to zero when t < 3 and are always less than 10 −5 during the whole evolution period of 5[K] −1 .Note that EE displays a jump at t = 0 because we introduce a magnon-like excitation at the boundary.Then, the EE gradually increases until it saturates to a plateau value.

SUPPLEMENTARY NOTE 3: PARTON MEAN-FIELD THEORY
In this section, we provide a detailed description of the parton mean-field theory defined in Eq. ( 4) of the main text.We summarize the main results as follows: We decouple the t-K-Γ model in the FM phase into separate hole and spin parts, H ≈ H h + H s , using a mean-field theory.The hole part H h can be effectively described by H h in Eq. (22).For the spin part, the groundstate energy corrected by the zero-point quantum fluctuation is E s in Eq. (34), which leads to a quantum order-by-disorder phenomenon for Γ = 0.The corresponding magnon excitations can be effectively described by a spin-wave theory with Hamiltonian H s ≈ H (2)  s + δH (2)  s , as defined in Eqs. ( 31) and (36), respectively.As mentioned in the main text, we represent fermionic electron operators in terms of fermionic holons h j and bosonic spinons b j,s : This representation, which is particularly suitable for studying symmetry-breaking states, enlarges the local Hilbert space.In order to restore the original Hilbert space, one needs to impose the local constraints, h † j h j = δ and s b † j,s b j,s = 1 − δ, for every lattice site j.By means of Eq. ( 15), the kinetic terms can be rewritten as We introduce a vector of spinons, b † j = b † j,↑ , b † j,↓ , and express the bilinear spin-spin interactions in a compact form as Then, we can obtain the spin part Hamiltonian on the NN bonds ⟨ jk⟩ ∈ γ, as Assuming a translational invariant ansatz, the Hamiltonians for the hole and spin parts can be further decoupled into several quadratic terms by introducing mean-field parameters on the NN bonds ⟨ jk⟩ ∈ γ, which should be determined self-consistently from Here P γ represents the holon mobility.And generally, F γ and A γ represent the ferromagnetic and antiferromagnetic short-range correlations, respectively.After implementing the following decoupling scheme as the mean-field Hamiltonian for the hole part reads where µ is the Lagrange multiplier to tune the filling number of holons such that ⟨h † j h j ⟩ = δ.Since P γ ∼ δ, the terms in the last line in Eq. ( 22) can be neglected in the small doping limit.
The four-boson terms Ĥγ jk in the spin part Hamiltonian need to be further decoupled in order to achieve a quadratic form.For instance, we employ the following mean-field treatment for Ĥx jk , which enables us to decouple the Kitaev-type interactions as follows: and the anisotropic spin-spin interaction as Similar decomposition schemes are also applied for Ĥy jk and Ĥz jk .For convenience, we further define on a NN bond ⟨ jk⟩ ∈ γ.Since the system exhibits translational invariance, we can perform a Fourier transform on the bosonic spinons and introduce a spinon vector at momentum k as where the additional index l = 0, 1 denotes the two sublattices of a honeycomb lattice.Using D γ , ψ † k , and a Lagrange multiplier λ to tune the condition such that n b ≡ ⟨b † j,s b j,s ⟩ = 1 − δ, we arrive at the mean-field Hamiltonian for the spinon part where two 4 × 4 matrices T γ k and D γ k correspond to bosonic hopping and pairing terms.The large-N limit of H s is taken with a parameter [72] κ ≡ n b /N being fixed.For a small value of κ, the Schwinger-boson mean-field theory may give rise to a gapful QSL, while a magnetically ordered state with bosonic spinon condensation appears at large κ.In our study, we focus on its realistic limit with κ = n b = 1−δ, i.e., N = 1.Here the explicit forms of the terms in Eq. ( 25) are The quadratic Hamiltonians H h and H s can be diagonalized by eigenvalue decomposition and bosonic Bogolyubov transformation, respectively.Here we denote the dispersion relations for the Schwinger bosons (spinons) as ϵ k,n with n = 1, 2, 3, 4.This enables us to numerically obtain the self-consistent solution for Eqs. ( 19)-( 21) through an iterative scheme.For instance, single occupancy leads to where ρ denote the condensate density, and 4 × 4 matrices V 11 k , V 12 k , V 21 k , and V 22 k form the symplectic matrix M k diagonalizing the Hamiltonian in Eq. ( 25), as . When ρ > 0 is a finite value, the Schwinger bosons condense at specific momenta This leads to a long-range magnetic ordered state.For instance, a ferromagnetic order along the z direction means ⟨b † j,↑ ⟩ = ⟨b j,↑ ⟩ = √ ρ = √ κ and similarly, one along the x direction When ρ = 0, all eigenmodes of Eq. ( 25) are gapped and the corresponding ground state is a disordered quantum spin liquid state.When A γ = 0 and thereby V 12 k = 0, one must have ρ = κ = 1 − δ corresponding to a ferromagnetic ordered ground state.
When carrying out the calculations, different randomly generated mean-field parameters are selected to initialize the iterations.However, after careful investigations, we conclude that, in the parameter regime of FM Kitaev coupling and subdominant Γ, there is no gapped QSL solution.It is well-known from the exact solution that the Kitaev honeycomb model exhibits a gapless Z 2 QSL.Hence, it is not surprising that the Schwinger-boson mean-field theory fails to describe the QSL phase for a FM Kitaev honeycomb model and only predicts a ferromagnetic ordered state, which is equivalent to a spin-wave theory.In the following, we would like to adopt a ferromagnetic order ansatz, which allows us to solve the mean-field theory using spin-wave theory in a more straightforward manner. [001] [100] [010]

B. Spin-wave theory for the spin part Hamiltonian
Based on our DMRG results for t > K and the findings of Schwinger boson mean-field theory, we would like to focus on the FM ordered phases.Instead of the Schwinger boson representation, for simplicity, here we use Holstein-Primakoff (HP) representation.For a semi-classical ferromagnetic order along the (sin ϕ cos θ, sin ϕ sin θ, cos ϕ) direction, the "spin operators" (note that they are not the physical S = 1/2 spins) now can be expressed as Ŝ y j = sin θ sin ϕ(κ − 2a † j a j ) + sin θ cos ϕ κ − a † j a j (a † j + a j ) where a † j and a j are the creation and annihilation operators for the HP bosons (magnons), respectively.The Schwinger bosons are closely related to the HP bosons via their single occupancy constraints.When θ = ϕ = 0, it is easy to verify that by eliminating the b ↑ boson using the constraint, the correspondence is [1] b ↓ ↔ a and b ↑ ↔ κ − a † a.
For the generic case, we can obtain that b ↑ → e −iθ/2 cos(ϕ/2) κ − a † a − e −iθ/2 sin(ϕ/2)a, b ↓ → e iθ/2 sin(ϕ/2) κ − a † a + e iθ/2 cos(ϕ/2)a.Substituting Eq. ( 27) into the spin part Hamiltonian, we can obtain a spin-wave mean-field theory up to the fourth-order: Before discussing each term in the above Hamiltonian, it is worth emphasizing that the self-consistent conditions become Eqs.( 19) and ( 29) which are numerically solved in an iterative manner.When carrying out the calculations, different randomly generated mean-field parameters are selected to initialize the iterations.
Here H (0) s is a constant term, namely, the semi-classical ground-state energy, where ℜ[•] denotes taking the real part, N c is the number of total unit cell, and Note that D ≡ D x = D y = D z is ensured by the C 6 lattice rotational symmetry.The first-order term, H (1)  s , is linear in the HP bosons: Note that both H (0) s and H (1)  s are independent of Kitaev coupling K, and thus, for Γ = 0, there are no constrains on the values of ϕ and θ.This indicates an emergent O(3) manifold of degenerate ground states for the pure Kitaev honeycomb model in the semiclassical limit.With Γ > 0, a stable solution of ferromagnetic order requires H (1)  s to vanish, which, in turn, leads to a magnetic order along the [111] direction with ϕ = arccos(±1/ √ 3) and θ = ±π/4.H (2)  s is the quadratic term in the Hamiltonian.In the basis of where T k and D k are 2 × 2 matrices: We emphasize that when Γ > 0, the spin polarization direction must be fixed as ϕ = arccos(±1/ √ 3) and θ = ±π/4.It is worth noting that due to the lattice and time-reversal symmetries, in practice, we have P x = ℜ[P x ] = P y = P z .Therefore, we verify from the form of T k that the hole kinetic terms can play a role similar to the ferromagnetic Heisenberg interactions in the spin-wave theory with J ≡ −⟨P x ⟩t/κ.

C. Order-by-disorder
Note that the constant term in Eq. ( 31) arises from the zero-point energy contribution of the bosons, i.e., appears because we express this quadratic Hamiltonian in a bosonic BdG form.Therefore, the ground state energy of H s up to the second order corrections reads where ω k,n are the magnon dispersions.Note that ω k,n still depends on θ and ϕ even when Γ = 0, indicating that the O(3) degeneracy at the Kitaev point is lifted by quantum fluctuations.We find that a [001] (or equivalently [100] and [010]) magnetization is more energetically favorable than other directions, see Fig. 11, indicating the presence of a quantum orderby-disorder mechanism.D. Flat magnon band in H (2)   s When we diagonalize H (2)  s , we find that the lower magnon band, ω k,0 , is a nearly flat band with almost zero energy, see Figs. 12, regardless of the spin polarization directions.This flat band arises because the semiclassical ground state of the ferromagnetic Kitaev honeycomb model is not necessarily a ferromagnetic order but instead a classical spin liquid.In order words, one can construct numerous zero-energy excitations on top of a ferromagnetically ordered ground state in the semiclassical Kitaev honeycomb model, as demonstrated in Fig. 13.Though we use the [100] and [xy0] order as examples to illustrate the zero-energy excitations, the same mechanism is applicable for arbitrary [xyz] order.Indeed, one can construct the same type of excitations on top of a [xyz] order by fixing the z components of all spins uniform and unchanged.Therefore, one can always observe flat magnon bands for arbitrary [xyz] ordered ground state in the pure Kitaev honeycomb model.
The flat band we discussed above can gain dispersion in two scenarios: (i) a finite positive Γ which also forces the polarization direction to be [111], and (ii) a large t with a finite doping level δ (i.e., a nonzero P γ ).Note that the lower band of H (2)  s exhibits gapless Goldstone-like modes even when Γ > 0, as shown in Fig. 14(a), although the system does not have any continuous symmetries to be spontaneously broken.This phenomenon has also been reported in Refs.[40,60].The zero-energy modes disappear in the presence of finite doping δ, as shown in Fig. 14(b).
Note that in this context, we have neglected the fourth-order contributions from the kinetic terms since P γ ∼ δ is relatively small in the small doping limit.Furthermore, the cubic terms vanish for the [001] order.We also find that the correction of cubic terms are not important for the [111] order, also see Ref. [75].
By introducing several Hartree-Fock parameters: the terms in H (4)  s on a α-type bond ⟨(r, 0), (r ′ , 1)⟩ are decomposed (r is the unit-cell index) as Then the quartic interaction terms can be replaced with where δE (4)  s and δH (2)  s are the Hartree-Fock corrections to the ground-state energy and magnon spectrum, respectively.
Substituting the above decomposition into the quartic terms in H (4)  s leads to Here the explicit form of δT k is Similarly, And the fourth-order correction to the ground-state energy is In the end, we would like to comment that the inclusion of δH (2)  s leads to several important effects in the spin-wave theory.In Fig 15, we show the magnon bands corrected by the fourth-order term δH (2)  s .We find that δH (2)  s significantly changes the dispersions of both the lower and upper magnon bands.In particular, the lower band gains much more dispersion, as shown in Fig. 15(a) and (b).Furthermore, the energy gap between the lower and upper magnon bands becomes closed due to the fourth-order corrections when Γ vanishes, see Fig. 15(b).It indicates that the topological properties of the magnon bands are also changed by δH (2)  s .Finally, a finite Γ can reopen the energy gap, as shown in Fig. 15(c).

SUPPLEMENTARY NOTE 4: ZIGZAG ORDER
Motivated by the experimental Kitaev materials, we discuss the effect of hole doping on the zigzag ordered state.Previous studies [29,46,48] report that the K-Γ model alone is insufficient to support a zigzag order, even with a dominant antiferromagnetic Γ < 0. Our simulations are consistent with these findings, as we observe something like spiral and/or incommensurate orders rather than zigzag order up to Γ = −5K, see Fig. 16(a).On the other hand, the t-K-Γ model at a doping level δ ≈ 0.01 can stabilize the zigzag ordered phase with a moderate t ≈ 0.2K with Γ = −5K, as shown in Fig. 16(b).This is attributed to that the hole kinetic terms act as a ferromagnetic Heisenberg exchange effectively, J ∼ |t|δ, which can prompt a zigzag order in the K-Γ model [66].
We find that a relatively large value of |Γ| > K is required to achieve a zigzag order in the t-K-Γ model, which might be due to the usage of narrow cylindrical geometries in the DMRG simulations [41].We anticipate that this value will significantly decrease in the two-dimensional limit.Alternatively, we introduce a small third NN antiferromagnetic Heisenberg exchange J 3 < 0 on top of the t-K (Γ = 0) model, which has shown to be sufficient to drive the system to a zigzag ordered phase [32,34,47], even on narrow cylinders.We find that in this scenario, a small t can also enhance the stability of the zigzag phase.For instance, when J 3 = −0.07K, a doping level of δ ≈ 0.01 and a hopping strength of t = 0.2K can approximately double the magnitude of zigzag order compared to the case at t → 0. Overall, our results suggest that a zigzag order can be more easily stabilized in the K-Γ (K-J 3 ) model by supplementing it with a moderate hole kinetic term t as well as small hole doping.Nevertheless, it is worth noting that if the hopping strength t dominates over the antiferromagnetic interactions of Γ < 0 (J 3 < 0), the system still undergoes a ferromagnetic ordering transition.Here we clarify the absence of a free spin moment in the gapless phase even in the limit of a static hole, i.e. hopping t = 0.As discussed in Ref. [49] of the main text, the KHM with a vacancy is still exactly solvable with the help of the Majorana representation.Denoting three Pauli matrices as σ γ (γ = x, y, z), we introduce σ γ j = ic 0 i c γ i with c 0 i and c γ i being the matter and gauge Majorana fermions, respectively.The doped hole introduces a static vacancy to the KHM.Without loss of generality, we assume that this vacancy is located at the site i 0 of A sublattice.Then, there are three unpaired gauge Majoranas (namely, Majorana zero modes), c x j , c y k , and c z l , where j, k, and l are three NN sites of i 0 .And there must be one Majorana zero mode γ 0 as a superposition of matter Majorana fermions.
After a Gutzwiller projection those four Majorana zero modes can lead to two-fold ground-state degeneracy.Moreover, in the gapped A phase [14], γ 0 is spatially localized around the static vacancy, as shown in Fig. 17   A phase, those four modes can effectively form a free and local spin-1/2 moment, leading to a finite magnetization of a single free spin-1/2 in the ground state manifold.However, in the gapless B phase, γ 0 is no longer a local zero mode, see Fig. 17 (b).Therefore, although the system exhibits an odd number of spins and the ground states must be two-fold degenerate, there is not necessarily a free spin-1/2 moment.In stark contrast, this effective spin-1/2 degree of freedom, which exhibits zero-energy excitation, strongly interacts with the rest of the system.In this case, even though the ground state has an odd number of spins, it does not give rise to a magnetization of a free spin-1/2 moment.

FIG. 1 .
FIG. 1. Magnetization in the hole-doped Kitaev spin liquid (KSL).(a) Illustration of the t-K model with three types of bonds, γ = x, y, z.Sites j = 1, ..., 6 label the convention for plaquette operators W p .The system sustains the site-diluted KSL ground state with fractionalized spins (W p ≈ 1) only for a slow hole (red ball).(b) For t ⪆ K doping depletes plaquette operators W p while simultaneously polarizing the spins, which leads to FM order along the [001] direction.(c,d) The ground-state phase diagram for the one-hole doped t-K model (K = 1) on a YC4 cylinder with length L x = 12 (corresponding to a doping level of δ ≈ 0.01).The evolutions of (c) the magnetic moment M z and (d) the average flux W p indicate a phase transition from a sitediluted KSL (blue regime in insets) to a FM state (red regime) around t * ≈ 0.28.

FIG. 2 .
FIG. 2. The magnetization as a function of doping level δ for variable values of hopping strength on YC4 cylinders with length L x = 12 and K = 1.

FIG. 3 .
FIG. 3. Effect of ferromagnetic off-diagonal symmetric exchange Γ.(a) The polar angle ϕ ≡ cos −1 (M z /|M|) representing the magnetization direction as a function of Γ, where ϕ = 0 indicates the [001] direction and the black dashed line the [111] direction.Inset: a finite FM Γ > 0 can induce a rotation of the spin polarization direction from [001] (transparent purple arrows) to [111] (solid red arrows) by an angle ϕ.(b) The amplitude of magnetization as a function of Γ for the one-hole doped model with t = 20K on a YC4 cylinder with length L x = 12, corresponding to doping level δ ≈ 0.01.Inset: The first-order derivatives of ground-state energy versus Γ.

FIG. 4 .
FIG. 4. Distinct boundary magnon excitations for the [001] and [111] polarized FM orders.(a) Magnon bands on a cylinder for a [001] ordered state within nonlinear SWT.(b, c) The real-time dynamics of (b) spin current j s (t, x) and (c) charge current j c (x, t) for a [001] ordered state.These results are obtained by DMRG using one-hole doped YC3 cylinders (δ ≈ 0.013).For (a, b, c), we set the offdiagonal exchange Γ = 0 and hopping strength t = 10K.The selfconsistent solution with doping level δ = 0.01 gives kinetic hole energy ⟨ Pjk ⟩ ≈ −0.93 × 10 −2 and FM spin correlation ⟨ F jk ⟩ ≈ 0.79 on the x− (y-) type bonds as well as ⟨ F jk ⟩ ≈ 0.82 on the z-type bonds.(d) Magnon bands within nonlinear SWT on a cylinder for the [111] ordered state.The chiral edge states are marked in green.(e, f) The real-time dynamics of (e) spin current j s (t, x) and (f) charge current j c (x, t) for a [111] ordered states within DMRG.For (d, e, f), we set Γ = 0.05K and t = 10K.The self-consistent solution with δ = 0.01 gives kinetic hole energy ⟨ Pjk ⟩ ≈ −0.93 × 10 −2 and FM spin correlation ⟨ F jk ⟩ ≈ 0.85 on all three type bonds.

SUPPLEMENTARY NOTE 1 :FIG. 6 .
FIG.6.Schematics of the labeling scheme for a honeycomb lattice as a cylindrical geometry with primitive vectors ⃗ x and ⃗ y.The blue, green, and red lines represent the x-, y-, and z-type bonds for a Kitaev honeycomb model, respectively.The periodic boundary condition has been imposed along the y direction, as indicated by the red dashed lines.

( 28 )
We have b † ↑ b ↑ + b † ↓ b ↓ = κ and hence κ = n b = 1 − δ is the density of spinons.Naturally, the short-range ferromagnetic correlation, up to the second-order of HP bosons, becomes

FIG. 13 .
FIG. 13.(a) In a [100] ordered ground state, where the effect of z-type (red dashed line) and y-type (green dashed line) bonds vanishes effectively, a zero-energy excitation can be constructed by flipping two spins (highlighted by the red shadow) connected by one x-type bonds (blue solid lines) into the [ 100] ( 1 = −1) direction.(b) Arbitrary numbers of zero-energy excitations can be created, leading to a flat magnon band for the [100] polarization.(c) In a [xy0] ordered ground state, the effect of z-type bonds (red dashed line) vanishes effectively and the honeycomb lattice is reduced into several independent chains formed by x-and y-type bonds (blue and green solid lines).The mechanism described in (a) and (b) no longer applies to this case.Instead, a zero-energy excitation can be created by rotating all the spins on a single chain (highlighted by the yellow shading) by an arbitrary angle around the z-axis.(d) All the spins on each decoupled chain can be rotated freely around the z-axis without costing energies, leading to numerous zero-energy excitations and thereby a flat magnon bond for a [xy0] polarization.Note that the mechanism described by (c) and (d) is also applicable to the [100] ordered state.

3 FIG. 16 .
FIG. 16.Zigzag order induced by mobile holes.Real-space magnetization ⟨S z j ⟩ (⟨S x j ⟩ and ⟨S y j ⟩ are similar) profiles for (a, b) t-K-Γ model on YC3 cylinders and (c,d) t-K-J 3 model (Γ = 0) on YC4 cylinders.Here J 3 is the 3rd NN antiferromagnetic Heisenberg exchange, K = 1, and the doping level is δ ≈ 0.015.Positive and negative magnetic moments are shown by red and blue dots, respectively, whose radius scales with the magnitude of magnetization |⟨S z j ⟩|, such that |⟨S z j ⟩| ≈ 0.06 in the two middle columns in (b).

FIG. 17 .
FIG. 17.The single-particle wave function of the Majorana zero mode γ 0 in the matter sector on a L y = 6 and L x = 12 honeycomb lattice.Note that we plot the honeycomb lattice as a 2L y × L x square lattice, where factor 2 accounts for the A and B sublattices.(a) The gapped A phase with K z = 3K x = 3K y .(b) The gapless B phase with K z = K x = K y .