Single-hole spectra of Kitaev spin liquids: from dynamical Nagaoka ferromagnetism to spin-hole fractionalization

The dynamical response of a quantum spin liquid upon injecting a hole is a pertinent open question. In experiments, the hole spectral function, measured momentum-resolved in angle-resolved photoemission spectroscopy (ARPES) or locally in scanning tunneling microscopy (STM), can be used to identify spin liquid materials. In this study, we employ tensor network methods to simulate the time evolution of a single hole doped into the Kitaev spin-liquid ground state. Focusing on the gapped spin liquid phase, we reveal two fundamentally different scenarios. For ferromagnetic spin couplings, the spin liquid is highly susceptible to hole doping: a Nagaoka ferromagnet forms dynamically around the doped hole, even at weak coupling. By contrast, in the case of antiferromagnetic spin couplings, the hole spectrum demonstrates an intricate interplay between charge, spin, and flux degrees of freedom, best described by a parton mean-field ansatz of fractionalized holons and spinons. Moreover, we find a good agreement of our numerical results to the analytically solvable case of slow holes. Our results demonstrate that dynamical hole spectral functions provide rich information on the structure of fractionalized quantum spin liquids.

Quantum spin liquids have long captivated interest as strongly-correlated systems due to their intriguing properties, such as long-range topological entanglement and fractionalized excitations [1][2][3][4][5].Theoretical models and classifications have been developed to understand these exotic states [6,7].The Kitaev spin liquid provides a unique route to access spin liquids because of its exact solvability [8].However, more relevant descriptions of materials include additional spin interactions spoiling exact solutions [9][10][11][12].Thus, mean-field analyses or numerical studies are required to explore their effects on the stability of the quantum spin liquid phase or the formation of competing magnetic orders [13][14][15][16].
Doping holes into the Kitaev spin liquid introduces an additional layer of complexity to the situation, enabling further fascinating phenomena.The model can still be solved exactly for slow holes by associating flux, fermion, and plaquette quantum numbers to the holes [17,18].Moreover, previous studies analyzed the influence of a finite hole density with additional spin terms in the form of Heisenberg interactions within mean-field theories [19][20][21][22][23]. Remarkably, they found several superconducting states, including possible topological superconductivity.
Another pending problem is the dynamics of holes in the limit of sparse doping or even for a single hole inserted into the spin-liquid ground state.This question is of particular interest for experimental probes such as angleresolved photoemission spectroscopy (ARPES) and scanning tunneling microscopy (STM), which directly measure the hole dynamics.Multiple plausible scenarios arise for the singlehole response.One possibility involves the fractionalization of the hole into holon and spinon quasiparticles.If these quasiparticles are deconfined, the hole spectral function will exhibit broad features and the spinons will determine the shape at the lowest energies at strong coupling.In a one-dimensional system, this scenario occurs naturally [24][25][26], but is also relevant for spectral functions in higher dimensions [27,28].This was directly confirmed numerically for the kagome spin liquid [29] and the chiral spin liquid on the triangular lattice [30].On the other hand, for confined spinons and holons, the spectrum shows a sharp, possibly renormalized, dispersion at low energies, as observed for the square lattice antiferromagnet [31][32][33][34][35][36][37][38][39][40].
Conversely, inserting a single hole can significantly impact the ground state of spin models.One famous example is the phenomenon of Nagaoka ferromagnetism [41,42], where the kinetic energy induced by the hole competes with the spin fluctuations.
For strong interactions or equivalently fast holes, the ground state may then favor ferromagnetic order.First directly observed in quantum dot experiments [43], this intriguing behavior has sparked recent explorations on frustrated lattices in the context of semiconductor heterostructures [44] and cold atoms [45][46][47].However, this phenomenon does not only affect single-hole ground states but can also emerge dynamically, as has been demonstrated for states at infinite temperatures [48,49].So far, dynamical Nagaoka ferromagnetism from quenching a ground state with a hole has remained unexplored.
Considering these possible scenarios, this work investigates the crucial question of what happens to a quantum spin liquid FIG. 1. Response of Kitaev spin liquids upon single-hole doping.When inserting a single hole into the Kitaev spin liquid, the hole either dynamically reorders the spin background into a Nagaoka ferromagnet for ferromagnetic (FM) spin couplings J > 0 or fractionalizes into spinon and holon quasiparticles for antiferromagnetic (AFM) spin couplings J < 0. arXiv:2309.15157v2[cond-mat.str-el]22 Mar 2024 state when a single hole is dynamically inserted.Specifically, we focus on the hole dynamics in the Kitaev spin liquid with both ferromagnetic and antiferromagnetic couplings; Fig. 1.As exact solvability is lost, we employ tensor network methods to simulate the time evolution after a hole quench and compute the hole spectral functions.In contrast to previous ARPES studies and high-temperature experiments of the Kitaev candidate materials α-RuCl 3 [50,51] and Na 2 IrO 3 [52][53][54], we directly analyze the fate of the Kitaev spin-liquid ground state upon sudden hole injection.Specifically, we show that for ferromagnetic Kitaev couplings, even slow-moving holes can dynamically polarize the ground state from a Kitaev spin liquid to a ferromagnetic state due to the dynamical Nagaoka effect.Intriguingly, for antiferromagnetic Kitaev couplings, the hole spectral function shows signatures of a spinon-holon fractionalization that a parton mean-field theory can partly explain.However, due to the complex interplay of spin, charge, and flux physics that renormalizes the energy scales of each other, the parton mean-field ansatz cannot capture all phenomena correctly.By contrast, in the limit of slow holes, we can directly confirm the contribution of fractionalized excitations [17].

Model
The exactly solvable Kitaev model [8] is a paradigmatic example of a quantum spin liquid hosting topological order and fractionalized excitations.Here, we summarize some of its most important properties.The Kitaev model is given by where the first, second, and third terms act on x-, y-, and z-bonds of a honeycomb lattice, respectively, see Fig. 1.
Throughout this work we will fix The system undergoes a phase transition from a gapped Z 2 topological phase for |J z | > 2 to a gapless phase for Remarkably, these phases exist for both signs of the Kitaev couplings J. Here, we will examine both ferromagnetic couplings J > 0 and antiferromagnetic couplings J < 0 and focus on the gapped case.We find qualitatively similar behavior for the gapless phase; see "Methods" section for additional data.The Kitaev model hosts an extensive number of conserved quantities, described by the flux operator around a plaquette; see inset in Fig. 2 (c): For each plaquette W has eigenvalues ±1.The ground state is in the flux free sector, where ⟨W⟩ = +1 for all plaquettes.Accordingly, we say that a flux is introduced into the system if ⟨W⟩ = −1 for one of the plaquettes.
For the exact solution of Eq. ( 1), the spin operators are decomposed into matter Majorana fermions χ 0 i and bond Majorana fermions Then, the spins are written in terms of the Majorana fermions as Introducing directed bond operators û⟨i j⟩ a = iχ a i χ a j from the A to the B sublattice, we rewrite the Kitaev model The gauge operators û⟨i j⟩ a commute with each other and with the matter Majorana fermions.Thus, the Hamiltonian splits into separate sectors of the bond operators, which have eigenvalues u ⟨i j⟩ a = ±1.Moreover, the plaquette operator W can be expressed as a product of û⟨i j⟩ a around a hexagon as W = û⟨12⟩ x û⟨23⟩ y û⟨34⟩ z û⟨45⟩ x û⟨56⟩ y û⟨61⟩ z .
To introduce holes into the system, two different approaches have been discussed before.On the one hand, in Refs.[19][20][21][22]55] holes are inserted in a t − J model formalism, similarly to the description of cuprates from a spin Hamiltonian on a square lattice [56].On the other hand, Refs.[17,18] introduce holes as spin sites, where the interactions to all other sites are missing.While the second approach gives exact results for very slow holes, we will focus on the first one, which captures the experimentally relevant finite hopping regime.Later, we compare it to the slow-hole limit as well.Therefore, we add a kinetic term for nearestneighbor hole hopping to the Kitaev Hamiltonian: c † jσ (c jσ ) create (annihilate) fermions with spin σ and the Gutzwiller projector P GW excludes doubly occupied sites.The holes are related to the spin operators by , where σ a are the corresponding Pauli matrices.The spin anisotropy in the Kitaev model originates from the strong spin-orbit coupling of electrons to the dorbitals, for example, in Iridate materiales [9].Therefore, the spins above are not electron spins but rather hybridized pseudo-spins.Nevertheless, they can directly be related to the electron operators at low energy as in Eq. ( 5), which have an isotropic hopping; see derivation in the "Methods" section.
When including the term H t , the Hamiltonian is no longer exactly solvable and also yields different behavior for AFM couplings J < 0 and FM couplings J > 0. Note that we stick to the usual Kitaev conventions with a minus sign in front of the spin Hamiltonian.Mean-field approaches offer convenient means to obtain results in this limit [19][20][21][22]57]; however, they possess inherent limitations and dependence on physical assumptions.Therefore, we will resort to numerical techniques to solve the dynamics of Eq. ( 5).We note that exact diagonalization studies of a similar system have been performed earlier [55].These are restricted to small system sizes and can hence only access a limited number of momenta for the hole spectral function.Here we use tensor networks in the form of matrix product states (MPS), which allow us  6); obtained from matrix product states (MPS) at times τ = 0 [1/J] and τ = 10 [1/J] are shown.We identify the FM cluster size N C by connecting all spins pointing up (purple dots) or down (orange dots) starting from the hole site (green dots).(b) Extracted from 10,000 snapshots of the state at each time, a large ferromagnetic (FM) cluster of spins forms around the hole.Of the N = 160 sites in the system, the cluster has a size of N c ≈ 35 (40) sites for hopping constants t/J = 3 (6).The speed of the FM cluster expansions is set by the hopping constant t; see dashed lines fitted to the snapshot data.(c) While the FM cluster size increases, the flux operator average ⟨W⟩ decays with time.All data shown is for ferromagnetic couplings J > 0 and anisotropy J z /J = 2.5.
to study larger systems on a cylinder with length L x and circumference L y .Details on the numerical methods are provided in the "Methods" section.

Dynamical Nagaoka ferromagnetism
The presence of a ferromagnetic (FM) Kitaev coupling constant (J > 0) is expected to be realized in some prominent Kitaev materials, such as Na 2 IrO 3 [58] or α-RuCl 3 [59].This motivates us to investigate the behavior of a single hole that is inserted into the spin-liquid ground state with FM couplings.
We first look at snapshots of the state during the time evolution and analyze these to shed light on the underlying dynamics.With the perfect-sampling algorithm introduced in Ref. [60], we can extract typical product-state configurations of the MPS and the corresponding probabilities in the computational basis, which is either a spin-up, spin-down, or hole state for the t-J model.For this basis {|s⟩ ≡ |s 1 s 2 . ..⟩ ; s i ∈ (↑, ↓, •)}, the wave function is given by Then the probability for a snapshot (i.e., the Fock space configuration |s⟩) is p(s) = |c s | 2 .We can use this for sampling the expectation value of an operator O with a subset S of basis states, where O is diagonal in the corresponding basis, Note that particle number conservation enforces precisely one hole per snapshot.Examples of such snapshots are displayed in Fig. 2 (a).For the initial state at τ = 0 [1/J], some small FM clusters are found, which are formed only from the spin background of the Kitaev spin liquid due to quantum fluctuations.At late times, e.g., τ = 10 [1/J], we observe larger FM clusters around the site where the hole is located.The anisotropic spin interactions J z /J = 2.5 used in our simulations favor FM correlations in z-basis, which is used as the Fock basis to sample the snapshots.Hence, the snapshots directly depict the fluctuating FM order.
The system undergoes a complex evolution over time.We define the FM cluster size N c around a hole as the number of aligned spins that directly connect to the site where the hole is for a given snapshot.Note that for the initial state, the hole will always be at the origin, where it was inserted.However, after some time, it will have spread to different lattice sites.Therefore, we cannot capture the emergent cluster size with (local) expectation values but have to resort to the analysis of snapshots.The average of 10,000 snapshots reveals that the cluster size increases linearly in time; see Fig. 2 (b).The velocity of this increase depends on the hopping parameter t, as demonstrated by the presented fits (dashed lines).Notably, for fixed system size, the maximum size of the FM clusters is proportional to the ground state magnetization M 0 = i ⟨ψ 0,1h |σ z i |ψ 0,1h ⟩ in the presence of one hole, which also is approximately proportional to t [61].However, due to the energy constraints within unitary time evolution, the FM clusters cannot reach the ground-state magnetization value M 0 as the hole excites the system.Instead, we find N C ∼ M 0 /2, with M 0 determined numerically.We also find that increasing t or doping a finite but small hole density instead of a single hole leads to the same magnetization for larger system sizes.
The dynamical formation of the FM cluster can be understood from the same argument as in Nagaoka's theorem for ground states with a single hole.For a hole that hops on a bipartite lattice, the kinetic energy is minimized when the interference of different hopping paths is constructive.In particular, this is the case when all spins are aligned, and hence the hole motion does not change the spin pattern.The minimization of kinetic energy usually competes with antiferromagnetic spin interactions, which prefer antialigned spins.For example, for the square lattice AFM, a large hopping strength is required to establish Nagaoka ferromagnetism.There, numerical computations find a large FM polarization around the hole for t Nagaoka /J ≳ 30 [62].In contrast, for the Kitaev model, the local spin interaction is already ferromagnetic and frustration only arises between the different x−, y−, and z− spin directions.Thus, a small hole hopping t/J and anisotropy J z can easily favor the ferromagnetism through kinetic energy minimization.This is FIG. 3. Hole spectral function for ferromagnetic Kitaev couplings.The spectral function A(k, ω) for ferromagnetic J > 0, J z /J = 2.5 and t/J = 3.0 along the three distinct cuts in the Brillouin zone (see inset in the bottom row) shows significant resemblance with the spectrum of a free hole hopping on a honeycomb lattice (green lines).Both the energy dispersion (top row) and the spectral weight (bottom row) agree well with the free hole.
the reason for the formation of large FM clusters around the hole, already for the moderate hopping strength of t/J = 3(6) considered in Fig. 2.
The flux average ⟨W⟩ over all plaquettes of the system, Fig. 2 (c), quickly decreases with time which indicates the disappearance of the Kitaev spin liquid.This observation highlights the transformative nature of the system as it transitions from a spin liquid to a ferromagnetic state.However, finite flux ⟨W⟩ can still be present simultaneously with large FM clusters, which is a remnant of the Kitaev spin liquid.Even though the spins get dynamically polarized, information remains about the initial flux configuration.The same phenomenon occurs in the ground state of the onehole doped Kitaev model, where both the magnetization and the flux expectation value ⟨W⟩ can be non-zero simultaneously [61].
The hole spectral function, as measured by ARPES, provides energy-and momentum-resolved information about the response of the spin-liquid state.It gives insights into the fractionalization of quasiparticles in terms of spinons and holons [28][29][30]63].Concretely, we are interested in calculating where |ψ 0 ⟩ denotes the Kitaev spin-liquid ground state.For numerical details and additional data regarding convergence in bond dimension, we refer to the "Methods" section.We shift all energies by a constant µ that is computed from the energy difference between the ground state without a hole and with a single hole.Therefore, it corresponds to the Mott gap of the insulator, and all hole excitations have to be outside that gap.In our convention, we plot −ω + µ, which hence must be a positive energy.
In previous X-ray scattering studies [64], the bonddepending interactions could be directly observed by measuring different spin components separately.Similarly, a spin-resolved response of the hole spectrum is experimentally feasible [65].We derive the connection between the spectrum for all physical electrons in a d 5 orbital and the low-energy effective electrons as in Eq. ( 8) in the "Methods" section.However, in our numerical computations, we observe the same responses for spin up and down for all spectral functions below, corresponding to the symmetry of c † j↑ and c † j↓ in Eq. (5).Therefore, we cannot gain additional insights from looking at the spin-resolved spectra numerically.However, spin-asymmetry may be used in experiments to characterize perturbations to the considered model.
The resulting spectrum for J > 0 with a Kitaev anisotropy J z /J = 2.5 and a hopping constant of t/J = 3.0 is shown in Fig. 3.We observe that most parts of the spectral function resemble that of a free hole (green line), meaning that it follows the dispersion from only H t in Eq. ( 5).This is consistent with a large FM cluster formation, where the hole can hop ballistically without distorting the spin background.Notably, both the energy dispersion and distribution of spectral weight Z k , obtained from this simple model, exhibit good agreement with the full spectrum of the Kitaev t − J model.For the spectral weight, we fit a Gaussian to the low energy (−ω + µ < 3t) and high energy (−ω + µ > 3t) branch of the MPS spectrum and integrate over the energies.The direct comparison of the spectral weight shows that indeed there are two branches in the spectrum, as for the free hopping case, but at each momentum, most of the weight is found in only one of them.Additionally, a more careful look at the spectral weight reveals missing weight in the sum of low and high energy branches for the MPS data, even though the sum rule dωA(k, ω) = 1 is fulfilled, when numerically integrating over the whole spectrum.This indicates that the spectrum also has some contributions from a continuum, which is especially visible in the M 2 -M 3 cut of the spectrum.There, the broadening is much wider than the applied Gaussian broadening for the numerical algorithm.Hence, even though the free-hole picture accounts for most features in the spectrum, it does not fully describe all the details, where the continuum suggests further interactions between several individual constituents.Although a broad continuum is expected when the hole separates into fractional quasiparticles, these features here are too fragile to infer the spin-liquid nature of the ground state from the hole spectrum.By contrast, the hole spectrum directly uncovers the dynamically emerging ferromagnetism induced by dynamically doping the hole, which is an interesting phenomenon on its own.
To further investigate the connection between the flux decay, as shown in Fig. 2 (c), and the hole motion, we look at the correlation functions ⟨n r W p ⟩ for a fixed hole position r six sites away from the original site of the hole insertion at time τ 0 = 0 [1/J].Initially, the hole is positioned at a fixed site in the middle of the cylinder, resulting in an expectation value of zero for the hole everywhere else and hence leads to a vanishing correlator.As time progresses, the hole propagates towards site r along various paths, as illustrated in Fig. 4 (a).After a short time τ 1 ≈ 0.5 [1/J], the shortest path that is a direct connection from the initial sites reaches site r, leading to the reduction of ⟨W p ⟩ on plaquettes p along that particular path.Subsequently, at longer times τ 2 and τ 3 , additional paths become accessible, resulting in the decay of ⟨W p ⟩ on more plaquettes in a wider range around the sites involved.Since the flux W p decays in a large area between the hole and its initial position, we can disregard the possibility of a composite quasiparticle, such as bound flux-hole states, that could move around freely but restores the spin background after propagating further.
To summarize, for ferromagnetic Kitaev couplings J > 0, the hole does not separate into fractionalized quasiparticles.Instead, the spin liquid state is destroyed, paving the way for a dynamically emerging ferromagnetic state.In the thermodynamic limit, only an infinitely fast hole is expected to give rise to an extensive FM polarization of the system.Our analysis, however, shows that already for intermediate hopping strength significant FM polarization clouds can be observed.Above, we focus on the specific gapped case of J z /J = 2.5 and t/J = 3.0.The same phenomena also occur for other ratios of t/J; data are shown in the "Methods" section.
Antiferromagnetic Kitaev couplings: Fast holes Antiferromagnetic (AFM) Kitaev couplings J < 0, which are predicted to be relevant for a recently explored class of Kitaev materials [66], exhibit distinct behavior upon doping compared to ferromagnetic (FM) couplings, as demonstrated at the mean-field level [21,22].Thus, we also investigate the dynamics of a hole inserted into a Kitaev spin-liquid ground state with antiferromagnetic couplings.For very large values t/|J| → ∞, Nagaoka's theorem is also applicable in this case, and an FM state will arise.However, as for the square lattice Heisenberg AFM, the hopping constant t required for that may be very large because the gain in the kinetic energy of the hole has to be balanced with AFM spin fluctuations.Here, we focus on the strong coupling regime t > |J| but avoid the Nagaoka FM for all considered values of the hopping strength.
In contrast to the FM case, the spectrum for AFM couplings exhibits a strikingly different character, Fig. 5 (a).There is no clear signature of free hole hopping.Instead, the presence of several flat bands suggests the importance of Kitaev physics and the localization of excitations.At low energies (−ω+µ) < 1.5t, the spectrum reveals a rich structure of dispersionless as well as dispersive bands; Fig. 5 (b).
Analytically understanding the dynamics of AFM couplings proves to be challenging due to the intricate interplay between the hole, spin, and flux physics.A possible approach to handle these different degrees of freedom is by fractionalizing the hole into a spinon and holon.Spinons couple to the spin degrees of freedom and hence contribute to the low-energy physics, while holons, carrying the charge quantum numbers, exhibit faster motion and contribute to high-energy features.The spectral function arises from the convolution of these two contributions [27,28]: On a one-dimensional lattice, the factorization of the spectral functions comes intrinsically from the spin-charge separation by mapping the t-J Hamiltonian to independent spinon and holon parts in squeezed space [63].By contrast, in two dimensions, the decoupling of spinons and holons from the hopping term generally can only be applied on a meanfield level.This approach yielded promising results for the chiral spin liquid [30].There, the lowest edge of the spectrum directly agrees with the dispersion of the corresponding spinon mean-field theory.Applying the same ansatz to the Kitaev spin liquid, we first rewrite the holedoped Hamiltonian Eq. ( 5) in terms of holons and spinons, which interact with each other.Then, we perform a meanfield decoupling to extract the separate spectra of these fractionalized excitations; see "Methods" section for more details.
The convolution of spinons and holons features several flat bands similar to the MPS spectrum.These can be identified from the flat spinon bands corresponding to bond excitations of the χ x,y,z Majorana fermions.For the anisotropic case J z /J = 2.5, the z-bands are at higher energy than the xand y-bands.The holon contributes to the convolution with a renormalized free hopping, which has the most spectral weight at the two van Hove singularities.Therefore, the most prominent signatures in the combined spectrum are flat bands at energies of the spinon bands plus holons at the van Hove singularities.As a reference, we also put these energy lines into the low energy spectrum obtained from MPS; see dashed and dotted lines in Fig. 5 (b) for the x-/y-spinons and z-spinon, respectively.
Overall, the parton mean-field theory captures the correct energies at which responses in the low-energy part of the hole spectrum are expected but fails to predict the spectral weight distribution; see Fig. 11 in the "Methods" section.Therefore, the mean-field ansatz does not encompass all the relevant correlations and interplay between the charge and spin degrees of freedom.This shortcoming may be associated with the strong holon renormalization leading to effective hoppings t z eff ≈ 0.18 t along the z-bonds and t x,y eff ≈ 0.11 t along the xand y-bonds; see "Methods" section for more details.Hence, the energy scales of spinons and holons are mixed but may be separated for even more extreme ratios t/|J|, which we cannot reliably access numerically.
To further compare the AFM couplings with FM couplings, we define the spectral density by integrating the spectrum over momentum D(ω) = 1 L x L y k A(k, ω); see Fig. 5 (c).Each of the spectra is cut off at the edge of the free hole motion at 6 t up to Gaussian broadening.However, while for FM couplings, the free hole dominates in the spectrum directly as a clear branch, for AFM couplings there are also flat bands in the high-energy regime of the spectrum, suggesting interactions of the free hole part with the flat band spinons, similarly to the renormalized holon at low energies.Another striking difference is that for AFM couplings, a part of the weight is shifted towards lower energies, where there is almost no spectral weight in the FM case.The absence of weight for FM couplings is consistent with the free-hopping picture.Nevertheless, the total spectral density has some deviations from the one expected from only the free hole.Instead of expected sharp peaks at the van Hove points, −ω + µ = 2t and −ω + µ = 4t, the signal is much more broadened, possibly due to further interactions between the hole and other constituents.The low-energy features present for the AFM couplings indicate that here the free hole does not describe the spectrum.To further support the absence of a FM polarization cloud around the hole, we conduct the same snapshot analysis as for the FM couplings.As shown in Fig. 5 (d), the cluster size increases over time in the case of FM couplings.However, the cluster size remains constant at a low value for AFM couplings.
In the "Methods" section we present data for different hopping amplitudes t/|J|.The spectra look very similar for different hoppings and only differ slightly in the distribution of spectral weight.Similar to FM couplings, also for AFM couplings the overall energy scaling of all detected signatures is proportional to t rather than J.This suggests a hole-dependent nature of these features and not only spin physics.By contrast, for the chiral spin liquid on the triangular lattice, the main energy scaling was found to be ∝ J [30].These observations let us conclude that for AFM couplings, holes, spins, and fluxes mutually influence each other, leading to the intricate dynamical interplay between the (fractionalized) excitations.The main reason for that is the strong renormalization of the hopping strength t eff , which is no longer larger than the spin energy scales.Nevertheless, the hole spectrum shows signatures of flat spinon bands that are a direct consequence of the underlying Kitaev spin liquid and can therefore be used as a characteristic of this phase in ARPES experiments.
Antiferromagnetic Kitaev couplings: Slow-hole limit Previous works have focused on the slow-hole limit t ≪ |J| as well [17,18].For a stationary hole t = 0, one can think of it as a vacancy site, where holes are described within the spin Hamiltonian Eq. ( 1), i.e., the spin is not removed from the "hole site", but all interactions with neighboring spins are just turned off.We now want to investigate this limit to understand the limitations and merits compared to our numerical simulations.
We start with the description of a stationary hole t = 0 and solve the quadratic Hamiltonian Eq. ( 4) by doubling the number of Majorana fermions to obtain complex fermions: where χ 0 i are copies of the original Majorana fermions [18].From this, we obtain the doubled Hamiltonian, We introduce a hole at site j 0 by switching off all interactions at the corresponding bonds, Let |Ω⟩ and | Ω⟩ denote the ground states of the pure Kitaev model and the one with a stationary hole, respectively.We can compute the corresponding quasiparticle weight of the hole spectrum Z = | ⟨ Ω|Ω⟩ | 2 directly on the same cylinder geometry that we use for the MPS time evolution.To obtain the full spectrum for a stationary hole, we expand FIG. 6. Slow-hole approximation.The spectral density D(ω) = 1 Lx Ly k A(k, ω) obtained from MPS time evolution (purple) is compared with the slow-hole approximation (green) as described in the main text and sketched in the inset of the last panel.We fix antiferromagnetic couplings J < 0, J z /J = 2.5 and vary t/|J| from 0.0 (exact stationary hole limit) to 3.0 (fast holes).For t/|J| = 3.0, a different energy scale of the x-axis is used in order to show the whole spectrum.
the hole Hamiltonian in the eigenbasis at half-filling, The spectral function is momentum independent since the hole is localized.In general, the sum above runs over exponentially many states.However, by evaluating the sum rule dωA(ω) = 1, we estimate that including only oneparticle excitations on top of the ground state already account for most of the spectral weight [67], e.g., for J z /J = 2.5, we find dωA 1 (ω) = 0.997.Therefore, when comparing this exact result to the MPS spectral function, we find excellent agreement up to little numerical deviations for the immobile hole; see Fig. 6.This serves as a stringent benchmark for our numerical approach.
To include a finite but small hopping amplitude t for the holes, we focus on the gapped phase J z /J > 2. In that limit, the effective hole hopping was derived in Ref. [17].This derivation is strictly valid as long as the hole excitations do not couple to any bulk excitations of the undoped Kitaev model, i.e., for t/|J| ≪ (J/J z ) 4 (the flux gap).Then the holes are associated with flux, fermion, and plaquette quantum numbers.In our case of a single hole, the former two have to be trivial, while the plaquette quantum number p = 0, 1 implies a degeneracy for the ground state | Ω⟩.We find, however, that all the calculations presented below yield identical results for both values.In the anisotropic limit, the hopping amplitude gets renormalized along the different bonds as following t x = t/2, t y = t/2, and t z = t [17].This can be intuitively understood from the isolated dimer limit J z /J → ∞, where the ground state consists of dimers on the z-bonds Inserting either an up-or downhole on one lattice site will break up the corresponding dimer, e.g., |↑ •⟩ z .The hole can hop freely back and forth along the same dimer on the z-bond |↑ •⟩ z → |• ↑⟩ z .Hopping along the x-or y-bonds requires hopping to a neighboring dimer and therefore, effectively reduces the hopping amplitude by a factor of 1/2, because of the strong ferromagnetic energy penalty J z /J → ∞ in the first term after hopping.Away from the isolated dimer limit, the renormalization changes little as long as J z /J > 2 [17].
The doped Hamiltonian Eq. ( 13) is only exactly solvable for a specific hole configuration.To interpolate between the limit of a completely stationary hole and a free moving one in the translational invariant Kitaev ground state, a variational ansatz was put forward [18], H(ρ) = ρH d + (1 − ρ) H j 0 d ; sketched in the inset of Fig. 6.Here, we modify this ansatz and allow for the bond strength to vary over all lattice sites ρ → ρ({ j}).The interpretation follows straightforwardly from the stationary picture.As the hole moves through the system, at different times, different bonds will be modified depending on the probability of finding the hole at the corresponding site.According to the hole expectation value ⟨n j ⟩, we choose ρ j = 1 − ⟨n j ⟩.
To describe a local hole that spreads dynamically, we approximate the time evolution by step-wise adjusting ⟨n j (τ)⟩ and, hence, also ρ j (τ), where ⟨n j (τ)⟩ is determined by a free hopping hole with t x = t/2, t y = t/2, and t z = t according to the slow hole prediction.Therefore, we get a Trotterized time evolution for the response function with step-size δ = τ/N, where we introduce hole operators a ( †) j which locally switch off the couplings.Here, E 0 is the ground state energy of the undoped Kitaev model and for each H n (n = 1, . . ., N) we modify ρ as described before according to the corresponding hole expectation at time n • δ.Since all H n are quadratic fermionic Hamiltonians, we can diagonalize them By insertion of identities after each time step, the correlation function simplifies to a product of overlaps, In practice, again we do not need to sum over all eigenstates |α (n) i n ⟩: Because the Hamiltonians change only slightly from one time step to the other, it is sufficient to take the overlap with only the ground states |Ω (n) ⟩.We find that indeed the overlaps are very close to one.Thus, we can simplify the time-dependent correlations to The spectral function A(k, ω) is then computed by spatial and temporal Fourier transformations of the time-dependent correlation functions, similar to the evaluation of the data from the MPS time evolution.
We compare the slow-hole approximations to the spectra obtained from the full MPS time evolution in Fig. 6.For t/|J| = 0.0, the analytic description of stationary holes is exact, and we find that the spectra indeed are almost identical; demonstrating that our numerical results are well converged as discussed above.When increasing the hopping strength, we see a remarkable resemblance between the slow-hole limit and the MPS data.Note that the formally required limit for the approximation to be valid is t/|J| ≪ (J/J z ) 4 ≈ 0.025.This ensures that the hole does not couple to the flux excitation of the bulk.Yet, our numerics agree well up to t/|J| ≈ 0.5.This indicates that even at higher energy, the hole-flux interactions are not relevant to the shape of the spectrum.Instead, the hole is best described by vacancy sites that spread slowly over the whole system and modify the Kitaev spin-liquid ground state only slightly.Upon increasing the hopping strength t/|J| ≳ 0.8, the general shape between both curves is still similar, but small deviations between peak positions occur.Furthermore, for t/|J| ≳ 1.0 additional peaks in the MPS spectrum appear.This suggests that further modes in the model with significant spectral weight are excited.Since this can include complex interactions between the hole and flux or matter excitations in the Kitaev model, the simple ansatz for the slow holes cannot capture these features anymore.Eventually, at even faster hopping t/|J| ≳ 3.0, the spectra look very different, meaning that the slow-hole picture breaks down, as expected.We find that the overall bandwidth is not given by the renormalized hoppings t/2 as predicted from the slow hole ansatz but by the full bandwidth ∼ t.Moreover, we see a shift of spectral weight towards low energies, where distinct peaks are visible, which transition into a broad continuum at high energies, and the parton ansatz becomes more reasonable.

Local scanning tunneling microscopy spectra
So far, we focused on the momentum-resolved and energyresolved hole spectrum as measured by ARPES experiments, FIG. 7. Local spectra.The local spectral function S (ω), Eq. ( 18), as measured by scanning tunneling microscopy (STM), is shown for ferromagnetic couplings J > 0 and antiferromagnetic couplings J < 0; both with J z /J = 2.5 and t/|J| = 3.0.The inset focuses on the low energy regime that is compared to the parton ansatz (green line).
Eq. ( 8) or their sum rules.Here, we will consider the spatially averaged local hole spectrum, that is measurable with scanning tunneling microscopy (STM): Although for lattices with a single site per unit cell S (ω) is equal to the spectral density D(ω) = 1 L x L y k A(k, ω), this is not the case for the honeycomb lattice, which has a two-sublattice structure; see methods for details.Recent experimental proposals suggest that tunneling through a Kitaev spin liquid layer into a metal beneath will lead to a distinct response of the underlying spin liquid [68][69][70].
Here, we consider a different scenario, where tunneling occurs directly in and out of the spin liquid layer by injecting mobile holes.
The local spectra for both FM and AFM couplings stretch over a broad continuum bounded by the free-hole bandwidth of 6 t; see Fig. 7.In contrast to the ARPES results, the local spectra S (ω) of the FM possesses additional structure at low energies.These additional peaks cannot be explained with the free hole alone but indicate the additional dressing by local dynamical spin correlations.For AFM couplings, the low-energy peak agrees reasonably with the parton mean-field ansatz derived in the "Methods" section; see the green line in the inset of Fig. 7.This supports the concept of fractionalized quasiparticles in this energy region of the spectrum and offers a different probe to detect signatures of the Kitaev spin liquid.

Discussion
Our results reveal an intriguing behavior of a mobile hole doped into the Kitaev spin-liquid ground state.We looked at the real-time dynamics of the hole spreading and studied the energy and momentum resolved hole spectra, as measured by ARPES as well as the local spectral function, as obtained from STM experiments.Our simulations were carried out on finite cylinder geometries using tensor network methods, however, we expect at least qualitatively similar behavior for the 2D limit.
The ARPES response is drastically different for ferromagnetic (FM) and antiferromagnetic (AFM) Kitaev couplings between neighboring spins.The FM case is characterized by the formation of clusters of aligned spins leading to a dynamical Nagaoka ferromagnet that allows for coherent hole hopping.Intriguingly, the STM spectra show coherent dynamical features even at low energies, which are not present in ARPES.Conversely, AFM couplings robustly show multiple dispersionless features at low energies.In the slow-hole limit, we find that the hole does not couple to any other excitations and is described as a vacancy that modifies the spin liquid.However, when increasing the hopping strength, flux and matter excitations of the Kitaev model become relevant and significantly alter the spectral function.Understanding these modifications is challenging since the exact solution of the Kitaev model is lost and the hole interacts non-trivially with the spin background.A parton mean-field ansatz gives some insights into the spectral responses for AFM couplings but fails to reproduce the proper spectral weight along the different momentum cuts in the Brillouin zone.Thus, interactions beyond a parton mean-field theory are required in the considered parameter regime and may be investigated in more detail in future studies.
Our work shows that doping a spin-liquid ground state with a single hole can have several different outcomes.In one scenario, the hole can fractionalize into a slow spinon and a fast holon.In that case, the low-energy spectrum directly probes the spinon dynamics of spin liquids [27,28,30].In the other scenario, the hole significantly modifies the spinliquid ground state in its vicinity and dynamically forms a Nagaoka ferromagnet around the hole.Both scenarios are realized in the doped Kitaev model depending on the sign of spin interactions.This behavior should be contrasted to the spectral functions of Mott insulating states with magnetic order as, for example, described by the t − J model on the square lattice.On the one hand, the AFM hole spectrum exhibits a quasi-particle peak at low energies that can be effectively described by a bound state of a holon and a spinon [39].In contrast, for the Kitaev spin liquid, we assume the holon and spinon to be deconfined.On the other hand, a FM spin exchange J results in a FM ground state on the square lattice.The single-hole problem can be solved exactly since hole hopping leaves the spin background unmodified and describes the hopping of a free hole.This is in sharp contrast to our results for the Kitaev spin liquid, where the ground state does not have any order initially.Instead, the FM forms dynamically and, along with it, the effective free-holelike response in the spectral function.
It will be interesting to see in which category the response of other slightly doped spin liquids may fall.This also raises the question of the general stability of quantum spin liquids upon hole doping.Tensor network methods for ground states as well as dynamical response functions offer one possible route to study the structure of the phases that could be realized in the doped Kitaev model [61,71].Furthermore, they may shed light on which perturbations can stabilize topological superconductivity beyond the mean-field approach.

Physical spins and form factor
In the promising candidate materials exhibiting strong Kitaev interactions, such as 5d iridate compounds Na 2 IrO 3 and 4d α-RuCl 3 , the partially filled d 5 orbital is split into e g and t 2g (|xy⟩, |yz⟩, and |zx⟩) orbitals.Because of the strong spin-orbit coupling, the t 2g multiplet is further divided into a J eff = 3/2 quartet and a J eff = 1/2 Kramers doublet.The lowenergy physics is dominated by the J eff = 1/2 doublet with a reduced bandwidth.Here, we would like to demonstrate the relation between hybridized pseudo-spins and physical spins in realistic materials such as Na 2 IrO 3 .
The Kramers doublet, (c † ↑ , c † ↓ ), can be expressed in terms of t 2g orbitals as [72] where d † m,σ is an electron for orbital m ∈ (xy, yz, zx) and spin σ ∈ (↑, ↓) and the lattice site index has been omitted.We can introduce a vector of electrons and define a projector Then we obtain a compact form of (c † ↑ , c † ↓ ) = d † P. The physical spin operators at site j are defined as Now we would like to project the spin to the low-energy J eff = 1/2 bands.In order to do this, we first need to find the null space corresponding to the projector P, dubbed P ⊥ .Note that P ⊥ , whose explicit form is not important, can be obtained by extracting the eigenstates of zero modes of PP † .Indeed, those zero modes of PP † correspond to the basis of the high-energy J eff = 3/2 quartet, and we denote those modes as (c † 1 , ..., c † 4 ).Using P and P ⊥ , one can express the electrons as . Therefore, we can rewrite the physical spin operators in terms of c-fermions, and after projecting into the low-energy doublet sector, we find that The tensor of the dynamical hole spectral function is thus defined as FIG. 8. Convergence of the spectral function.The spectral function A(k, ω) for ferromagnetic couplings, J > 0, J z /J = 2.5 and t/J = 3.0 at the Γ point is converged for larger values of the maximal MPS bond dimension χ > 500, but still shows considerable deviations for χ = 300.
By viewing A(k, ω) as a 6 × 6 matrix, it is related to the spectral function for c-fermions by a unitary transformation as A(k, ω) = (P + P ⊥ ) * F (k, ω)(P + P ⊥ ) T , where, with H the Hamiltonian and E 0 the ground-state energy.We mainly focus on the excitations with energy being within E W , the bandwidth of the J eff = 1/2 band.Within the frequency regime of ω < E W (ℏ = 1), we approximately have Consequently, the orbital and spin resolved hole spectral function can be obtained by only computing the 2 × 2 spectral function for the pseudo spins F σ,σ ′ , as In the end, the energy-and momentum-resolved spectral function defined in the main text Eq. ( 8) is just the trace of the above matrix, as Hence, the only contributions come from terms that do not mix the different spins.Numerically, we find that F ↑,↑ (k, ω) = F ↓,↓ (k, ω) for all investigated parameters.Therefore, spin-resolved spectra provide only additional information when the material breaks spin-symmetry.

Details on the numerical methods
We use MPS methods for computing the time evolution of the state after hole injection.All simulations were performed with the Python library TeNPy [73].Numerical costs grow FIG. 9. Additional data for hole spectral functions with hopping amplitudes t/|J| = 6, and 12. (a), (b) For ferromagnetic Kitaev couplings (J > 0), the spectrum mainly consists of the response of a free hole (green lines).(c), (d) Antiferromagnetic Kitaev couplings (J < 0) show only slight changes at low energies for different values for t/|J|.Spectra are compared to the parton ansatz; dashed and dotted lines for spinon excitations along x-or y-bonds and z-bond, respectively.Note that all energy scales are in units of t, and we fix J z /J = 2.5.
linearly in L x but exponentially in L y .Thus, we restrict ourselves to cylinders with L y = 4 unit cells (i.e., eight sites) and L x = 20.To reach the large bond dimensions of up to χ = 1000 required to obtain converged results, we have utilized U(1) × Z 2 symmetries for particle number and spin parity conservation, respectively.
To study the hole dynamics, we first find the approximate ground state of the pure spin model Eq.(1) with DMRG.Several flux sectors are degenerate on the chosen cylindrical geometry with open boundary conditions along the xdirection.We favor the flux-free sector explicitly as a ground state in our simulation by adding plaquette terms ∝ − p W p .For a periodic or infinite system, it is known that the ground state has to be flux-free [74].
The time evolution was realized via an MPO representation of the time evolution operator, the W II operator from Ref. [75].To apply the operator to the state, we have first to use the more costly zip-up methods for a few time steps to avoid being stuck in a local minima before continuing with a variational truncation scheme with fixed MPS bond dimension χ.Depending on the hole hopping we choose step size δτ = 0.05/ t for fast holes t/|J| > 1.0 and δτ = 0.025/ |J| otherwise.
In this way, we compute a time evolution after injecting a hole into the ground state |ψ 0 ⟩ of the pure Kitaev model where j is a site in the middle of the cylinder.We check that the cylinder is long enough (L x = 20) such that excitations do FIG.10.Additional data for hole spectral functions with isotropic Kitaev couplings J z /J = 1.(a) For ferromagnetic Kitaev couplings (J > 0), the main contribution to the spectrum resembles the response of a free hole (green lines).(b) For antiferromagnetic Kitaev couplings (J < 0), the first flat features can be described by the parton ansatz (dashed green lines).We fix t/|J| = 3.0.
not reach the boundaries on the simulated time scales to limit the boundary effects.
To access the spectral function, we employ established MPS methods [76].
Initially, we obtain an MPS approximation for the ground state without hole |ψ 0 ⟩ through the DMRG algorithm.The calculation is performed on a finite L y × L x system, with L y = 4 and L x = 20.For that state, we compute the time-dependent correlation function Using translational invariance of the ground state, we only need to compute two time evolutions after hole insertion at one site in each of the A and B of the honeycomb lattice.For ARPES, the hole can be ejected by the photon on any of the sublattices.Therefore, the measured spectral weight is only periodic in an extended Brillouin zone [55,77] and the corresponding Fourier transformation of the annihilation operator is given by Here, the sum includes all lattice sites.This has to be contrasted to the usual definition of the sublattice Fourier transformation, which is for instance used to compute the dispersion relations for the parton mean-field ansatz below For Eq. ( 31), a unique inverse Fourier transformation exists.In contrast, for the operator in Eq. ( 30), we have to distinguish between k ∈ 1.BZ and k 1.BZ.These two cases are related to the sublattice Fourier transformation by: where g is the reciprocal lattice vector that brings k back to the first Brillouin zone and r AB is the vector connecting the two sites in the unit cell.We use these conventions to compute the spatial Fourier transformation of Eq. ( 29).The subtlety of contributions from both sublattices also leads to the difference between the local hole spectrum and the integrated spectral density S (ω) D(ω), see Eq. ( 18) in the main text.Furthermore, we can define a sublattice spectral function similar as in Eq. ( 8), but restrict the hole creation/ annihilation operator to one sublattice, Eq. (31).Although this spectrum can not be measured experimentally, it can be computed numerically and yields additional insights, as demonstrated for the Kitaev-Heisenberg model [55].In this sublattice spectral function, the lower branch of the freehole spectrum has finite spectral weight.By contrast, in the experimentally accessible hole spectrum Fig. 3, the lower branch remains dark due to cancellations from phase factors in Eq. (32).
To improve the finite time restrictions due to limited entanglement with a fixed bond dimension, we use linear prediction [78] combined with a Gaussian envelope after taking the spatial Fourier transform but before the temporal Fourier transform to obtain the spectral function, Eq. ( 8) in the main text.On the other hand, for the local spectrum, Eq. ( 18), we use i = j and sum over the contributions from the two sublattices without the spatial Fourier transformation but still keeping the linear prediction and Gaussian envelope.
To check convergence with bond dimension χ, we compute the time evolution with several values of χ ∈ {300, 500, 1000}.As seen in Fig. 8, some qualitative features already develop correctly for a small bond dimension χ = 300.However, the spectral weight still changes when increasing χ.From χ = 500 to χ = 1000, the spectral function remains unchanged, indicating convergence.Data shown in all other figures are obtained for bond dimensions χ = 500.

Additional data for spectral functions
We present additional spectra for different parameters.First, we investigate the spectra for different hopping parameters t when considering ferromagnetic (J > 0) and antiferromagnetic (J < 0) coupling constants.The resulting spectra are shown in Fig. 9.
For the ferromagnetic case, Fig. 9 (a)-(b), we observe only slight changes in the distribution of spectral weight with varying hopping parameters t.The general picture remains, further supporting that the spectrum is dominated by free hole hopping.All energy scales are given in units of t.Increased hopping leads to a higher contribution from the kinetic energy, which further stabilizes the ferromagnetic order and the associated free coherent hole motion.
Next, we consider the antiferromagnetic spin coupling, J < 0, and focus on the low-energy scales for different hopping parameters.We find that despite a small redistribution of spectral weight for different t, similar features are present throughout the spectra; see Fig. 9 (c)-(d).We compare the MPS data to the flat bands obtained from the parton construction.The energy difference between the anisotropic flat spinon bands scales with J z /J.Thus, the difference becomes smaller in units of t as we increase the hopping strength.Furthermore, for t/|J| = 12, the flat bands in Fig. 9 (d) are shifted towards slightly higher energies than predicted from the parton mean-field ansatz.This could indicate an additional interplay between spinons and holons.We note that even for these comparatively large values of the hopping, the effective hopping scale, which is renormalized down by a factor of five to ten, is still in the same order as the exchange.Thus, it is very difficult to numerically reach a true strong coupling limit in which the hole dynamics is effectively much faster than the spin dynamics.
Furthermore, we study the case of isotropic Kitaev J z /J = 1 in Fig. 10.This corresponds to the gapless phase of the Kitaev spin liquid [8].A larger MPS bond dimension is needed to faithfully represent even the ground state at half-filling.The subsequent time evolution is also more challenging, and we thus expect small deviations by further increasing the bond dimension beyond the accessible bond dimension of χ = 1000, which we used for these simulations.However, qualitatively, we can already see in Fig. 10 that the general behavior in the gapless phase is the same as in the gapped phase: the free-hole response dominates the spectrum for FM Kitaev couplings J > 0, while in the AFM case, J < 0 dispersionless features according to the parton picture are dominating.These features emerge from the flat spinon band together with the dominant holon spectral weight at the van Hove singularity.Here, we would expect a difference between the gapless and gapped phases.
While the gapped phase is characterized by flat spinon bands at different energies along the x-/y-bonds and z-bonds, they should have the same energy for isotropic couplings in the gapless phase.Unfortunately, we do not see this difference clearly in our numerical spectra, which could be due to the finite energy resolution or possible mixing of the spinons and holons when including further interactions.Moreover, similar to the gapped case, the flat spinons contribute dominantly to the spectrum, and the dispersive gapless spinon does not possess any significant matrix elements.We also do not see any other contributions in the MPS data that could be directly attributed to dispersive spinons.

Parton mean-field theory
When we add a hole hopping term H t to the Kitaev model, the whole Hamiltonian Eq. ( 5) is no longer exactly solvable.A general approach to describe spin liquids in the presence of doping is a parton mean-field ansatz.To compute the spectral function using this ansatz, we follow Ref. [19] and consider the splitting of a hole into a charge degree of freedom, the holon b † , and a spin degree of freedom, the spinon f † σ , Here, the holons fulfill the bosonic commutation relations [b in , b † jm ] = δ i j δ nm and the spinons obey fermionic anticommutation relations { f iσ , f † j σ} = δ i j δ σ σ.By introducing two holon species b † i1 , b † i2 we retain the SU(2) gauge redundancy.At each site we must impose the constraint The spinons can be related to the Majorana fermion description of the pure Kitaev model Eq. ( 4) by fixing a SU(2) gauge [14,79], Therefore, we can rewrite the Kitaev model in terms of these spinons: where the parameters J, ∆ are determined selfconsistently [79].
Next, we include the hole hopping term in our description.In our ansatz, we assume that the holon and the spinon are deconfined [27,28] and carry out a mean-field holons appear at different energy scales and the lower edge of the spectrum is given by the spinon dispersion [30,63].However, from Eq. ( 36), we see that the effective hopping constant is strongly renormalized.Concretely, for J z /J = 2.5, we find |W i j | ≈ 0.23 on x-and y-bonds and |W i j | ≈ 0.36 on z-bonds.Hence, the spin anisotropy directly results in an anisotropic hopping for the holon.With the prefactor 1/2 in Eq. ( 36), the hopping gets renormalized to quite small effective values.Thus, the spinons and holons no longer have separate energy scales, preventing a clear distinction between the two in the spectrum for reasonable parameters.
In Fig. 11, we compare the resulting spectra from the parton construction to the MPS data.The parton mean-field theory only gives predictions for the low-energy part of the spectrum.As can be seen in Fig. 11 (b) the convolution of holon and spinon according to Eq. (38) gives rise to several flat bands.This is expected from the spinon component, where flat bands arise from the bond Majorana fermions.Focusing on the anisotropic case J z /J = 2.5, we get distinct bands for the xand y-bonds (dashed lines) compared to the z-bonds (dotted lines).Importantly, the holon dispersion has a van Hove singularity, at which many states with the same energy will contribute to the spectral function convolution leading to an intense spectral peak.This gives rise to the flat band structure seen in the parton spectrum at energies ε h (k vH ) + ε x,y,z sp .At energies above the van Hove singularities of the holon, there are dispersing bands with very weak spectral weight, which originate from the convolution of the holon with the dispersive matter Majorana modes.
The MPS spectrum features similar flat bands as well; see Fig. 11 (a).The lower one is at the same energy as the parton construction and may be interpreted as the lowest spinon mode plus the van Hove singularity of the holon.However, the spectral weight distribution shows qualitatively different behavior, indicating that there are complex interactions between holon and spinon that must be treated beyond the simple mean-field level.At larger energies (−ω + µ) ∼ 1.0t, the MPS spectrum shows dispersive features as well, which have some resemblance of the parton ansatz, especially for the Γ-M momentum cut.However, the MPS spectra carry much more spectral weight at these energies than the parton predictions.

FIG. 2 .
FIG. 2. Dynamically emerging Nagaoka ferromagnetism in a doped Kitaev spin liquid with FM spin couplings.(a) Typical examples for snapshots of the Fock configurations occurring with large probability in the wave function; see Eq. (6); obtained from matrix product states (MPS) at times τ = 0 [1/J] and τ = 10 [1/J] are shown.We identify the FM cluster size N C by connecting all spins pointing up (purple dots) or down (orange dots) starting from the hole site (green dots).(b) Extracted from 10,000 snapshots of the state at each time, a large ferromagnetic (FM) cluster of spins forms around the hole.Of the N = 160 sites in the system, the cluster has a size of N c ≈ 35(40) sites for hopping constants t/J = 3(6).The speed of the FM cluster expansions is set by the hopping constant t; see dashed lines fitted to the snapshot data.(c) While the FM cluster size increases, the flux operator average ⟨W⟩ decays with time.All data shown is for ferromagnetic couplings J > 0 and anisotropy J z /J = 2.5.

FIG. 4 .
FIG. 4. Flux dynamics after hole insertion.(a) When hopping from the initial site to the response site r, the hole can move along different paths contributing to the dynamics at the corresponding times.(b)The correlations ⟨n r W p ⟩ between a hole at a fixed site r and all plaquettes p get destroyed along the paths as time evolves from τ 1 to τ 3 .Shown data is for ferromagnetic Kitaev couplings J > 0, J z /J = 2.5 and t/J = 3.0.

FIG. 5 .
FIG. 5. Hole spectral function for antiferromagnetic Kitaev couplings.(a) The spectral function A(k, ω) for antiferromagnetic couplings, J < 0, J z /J = 2.5 and t/|J| = 3.0 along the three distinct cuts in the Brillouin zone (see inset on the top right).(b) Zoom-in on the low-energy part of the spectrum shows several flat dispersions, which agree with the parton mean-field description (green dashed and dotted lines).(c) The spectral density D(ω) = 1 Lx Ly k A(k, ω), for antiferromagnetic (AFM J < 0) and ferromagnetic (FM J > 0) Kitaev couplings.(d) The FM cluster size N c increases for FM couplings but stays constant for AFM couplings, as further illustrated by the fits (dashed lines).