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

Abstract

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.

Introduction

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⁸. 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 angle-resolved photoemission spectroscopy (ARPES) and scanning tunneling microscopy (STM), which directly measure the hole dynamics. Multiple plausible scenarios arise for the single-hole 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²⁹ and the chiral spin liquid on the triangular lattice³⁰. 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⁴³, this intriguing behavior has sparked recent explorations on frustrated lattices in the context of semiconductor heterostructures⁴⁴ 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 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₃^50,51 and Na₂IrO₃^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¹⁷.

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.

Results

Model

The exactly solvable Kitaev model⁸ 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

$${H}_{K}=-{J}_{x}\mathop{\sum}\limits_{{\langle i,j\rangle }_{x}}{\sigma }_{i}^{x}{\sigma }_{j}^{x}-{J}_{y}\mathop{\sum}\limits_{{\langle i,j\rangle }_{y}}{\sigma }_{i}^{y}{\sigma }_{j}^{y}-{J}_{z}\mathop{\sum}\limits_{{\langle i,j\rangle }_{z}}{\sigma }_{i}^{z}{\sigma }_{j}^{z},$$

(1)

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 ∣J_x∣ = ∣J_y∣ ≡ ∣J∣ = 1. The system undergoes a phase transition from a gapped ${{\mathbb{Z}}}_{2}$ topological phase for ∣J_z∣ > 2 to a gapless phase for ∣J_z∣ < 2. 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. 2c:

$$W={\sigma }_{1}^{x}{\sigma }_{2}^{y}{\sigma }_{3}^{z}{\sigma }_{4}^{x}{\sigma }_{5}^{y}{\sigma }_{6}^{z}$$

(2)

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.

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.

For the exact solution of Eq. (1), the spin operators are decomposed into matter Majorana fermions ${\chi }_{i}^{0}$ and bond Majorana fermions ${\chi }_{i}^{a}$ [a ∈ (x, y, z)], such that $\{{\chi }_{i}^{\alpha },{\chi }_{j}^{{\alpha }^{{\prime} }}\}=2{\delta }_{ij}{\delta }_{\alpha {\alpha }^{{\prime} }}$ [α ∈ (0, x, y, z)]. Then, the spins are written in terms of the Majorana fermions as

$${\sigma }_{i}^{a}=i{\chi }_{i}^{0}{\chi }_{i}^{a},\quad a\in (x,y,z).$$

(3)

Introducing directed bond operators ${\hat{u}}_{{\langle ij\rangle }_{a}}=i{\chi }_{i}^{a}{\chi }_{j}^{a}$ from the A to the B sublattice, we rewrite the Kitaev model

$${H}_{K}=i\mathop{\sum}\limits_{a,{\langle ij\rangle }_{a}}{J}_{a}{\hat{u}}_{{\langle ij\rangle }_{a}}{\chi }_{i}^{0}{\chi }_{j}^{0}.$$

(4)

The gauge operators ${\hat{u}}_{{\langle ij\rangle }_{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}_{{\langle ij\rangle }_{a}}=\pm 1$. Moreover, the plaquette operator W can be expressed as a product of ${\hat{u}}_{{\langle ij\rangle }_{a}}$ around a hexagon as $W={\hat{u}}_{{\langle 12\rangle }_{x}}{\hat{u}}_{{\langle 23\rangle }_{y}}{\hat{u}}_{{\langle 34\rangle }_{z}}{\hat{u}}_{{\langle 45\rangle }_{x}}{\hat{u}}_{{\langle 56\rangle }_{y}}{\hat{u}}_{{\langle 61\rangle }_{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⁵⁶. 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 nearest-neighbor hole hopping to the Kitaev Hamiltonian:

$$\begin{array}{l}H={H}_{t}+{H}_{K}\\ \quad\,=-t\mathop{\sum}\limits_{\langle ij\rangle ,\sigma }{{{{\mathcal{P}}}}}_{{{{\rm{GW}}}}}\left({c}_{i\sigma }^{{\dagger} }{c}_{j\sigma }^{}+{{{\rm{h.c.}}}}\right){{{{\mathcal{P}}}}}_{{{{\rm{GW}}}}}-\mathop{\sum}\limits_{a,{\langle ij\rangle }_{a}}{J}_{a}{\sigma }_{i}^{a}{\sigma }_{j}^{a}.\end{array}$$

(5)

${c}_{j\sigma }^{{\dagger} }$ (${c}_{j\sigma }^{}$) create (annihilate) fermions with spin σ and the Gutzwiller projector ${{{{\mathcal{P}}}}}_{{{{\rm{GW}}}}}$ excludes doubly occupied sites. The holes are related to the spin operators by ${\sigma }_{i}^{a}=({c}_{i\uparrow }^{{\dagger} },{c}_{i\downarrow }^{{\dagger} }){\sigma }^{a}{({c}_{i\uparrow }^{},{c}_{i\downarrow }^{})}^{T}$, 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 d-orbitals, for example, in Iridate materiales⁹. 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⁵⁵. 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 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₂IrO₃⁵⁸ or α-RuCl₃⁵⁹. 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. ⁶⁰, 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 $\{\left\vert s\right\rangle \equiv \left\vert {s}_{1}{s}_{2}\ldots \,\right\rangle ;{s}_{i}\in (\uparrow ,\downarrow ,\circ )\}$, the wave function is given by

$$\left\vert \psi \right\rangle =\mathop{\sum}\limits_{s}{c}_{s}\left\vert s\right\rangle .$$

(6)

Then the probability for a snapshot (i.e., the Fock space configuration $\left\vert s\right\rangle$) is p(s) = ∣c_s∣². We can use this for sampling the expectation value of an operator O with a subset ${{{\mathcal{S}}}}$ of basis states, where O is diagonal in the corresponding basis,

$$\left\langle \psi \right\vert O\left\vert \psi \right\rangle \approx \mathop{\sum}\limits_{s\in {{{\mathcal{S}}}}}p(s)\left\langle s\right\vert O\left\vert s\right\rangle .$$

(7)

Note that particle number conservation enforces precisely one hole per snapshot. Examples of such snapshots are displayed in Fig. 2a. 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. 2b. 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}={\sum }_{i}\langle {\psi }_{0,{{\mbox{1h}}}}| {\sigma }_{i}^{z}| {\psi }_{0,{{\mbox{1h}}}}\rangle$ in the presence of one hole, which also is approximately proportional to t⁶¹. However, due to the energy constraints within unitary time evolution, the FM clusters cannot reach the ground-state magnetization value M₀ as the hole excites the system. Instead, we find N_C ~ M₀/2, with M₀ 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 anti-aligned 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⁶². 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 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. 2c, 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 one-hole doped Kitaev model, where both the magnetization and the flux expectation value 〈W〉 can be non-zero simultaneously⁶¹.

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

$$A({{{\bf{k}}}},\omega )=\mathop{\sum}\limits_{\sigma =\uparrow ,\downarrow }\int{{{\rm{d}}}}\tau {{{{\rm{e}}}}}^{i\tau \omega }\left\langle {\psi }_{0}\right\vert {c}_{{{{\bf{k}}}}\sigma }^{{\dagger} }(\tau ){c}_{{{{\bf{k}}}}\sigma }^{}(0)\left\vert {\psi }_{0}\right\rangle ,$$

(8)

where $\left\vert {\psi }_{0}\right\rangle$ 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⁶⁴, the bond-depending interactions could be directly observed by measuring different spin components separately. Similarly, a spin-resolved response of the hole spectrum is experimentally feasible⁶⁵. We derive the connection between the spectrum for all physical electrons in a d⁵ 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\uparrow }^{{\dagger} }$ and ${c}_{j\downarrow }^{{\dagger} }$ 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, ${Z}_{{{{\bf{k}}}}}^{(-\omega +\mu\, > \,3t)}+{Z}_{{{{\bf{k}}}}}^{(-\omega +\mu\, < \,3t)} < 1$, 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₂-M₃ 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.

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.

To further investigate the connection between the flux decay, as shown in Fig. 2c, and the hole motion, we look at the correlation functions 〈n_rW_p〉 for a fixed hole position r six sites away from the original site of the hole insertion at time τ₀ = 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. 4a. After a short time τ₁ ≈ 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 τ₂ and τ₃, 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.

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_rW_p〉 between a hole at a fixed site r and all plaquettes p get destroyed along the paths as time evolves from τ₁ to τ₃. Shown data is for ferromagnetic Kitaev couplings J > 0, J_z/J = 2.5 and t/J = 3.0.

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⁶⁶, 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. 5a. 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. 5b.

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(\omega )=\frac{1}{{L}_{x}{L}_{y}}{\sum }_{{{{\bf{k}}}}}A({{{\bf{k}}}},\omega )$, 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).

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:

$$A({{{\bf{k}}}},\omega )=\int{{{\rm{d}}}}\nu {{{\rm{d}}}}{{{\bf{q}}}}\,{A}_{{{{\rm{sp}}}}}({{{\bf{k}}}}-{{{\bf{q}}}},\omega -\nu ){A}_{{{{\rm{h}}}}}({{{\bf{q}}}},\nu ).$$

(9)

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⁶³. By contrast, in two dimensions, the decoupling of spinons and holons from the hopping term generally can only be applied on a mean-field level. This approach yielded promising results for the chiral spin liquid³⁰. 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 hole-doped Hamiltonian Eq. (5) in terms of holons and spinons, which interact with each other. Then, we perform a mean-field 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 x- and 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. 5b 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}_{{{{\rm{eff}}}}}^{z}\,\approx \,0.18\,t$ along the z-bonds and ${t}_{{{{\rm{eff}}}}}^{x,y}\,\approx \,0.11\,t$ along the x- and 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(\omega )=\frac{1}{{L}_{x}{L}_{y}}{\sum }_{{{{\bf{k}}}}}A({{{\bf{k}}}},\omega )$; see Fig. 5c. 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. 5d, 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³⁰. 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:

$${f}_{i\in A}=\frac{1}{2}\left({\chi }_{i}^{0}+i{\underline{\chi }}_{i}^{0}\right),\quad {f}_{i\in B}=\frac{i}{2}\left({\chi }_{i}^{0}+i{\underline{\chi }}_{i}^{0}\right),$$

(10)

where ${\underline{\chi }}_{i}^{0}$ are copies of the original Majorana fermions¹⁸. From this, we obtain the doubled Hamiltonian,

$${H}_{d}={H}_{K}+{\underline{H}}_{K}=i\mathop{\sum}\limits_{a,{\langle ij\rangle }_{a}}{J}_{a}{\hat{u}}_{{\langle ij\rangle }_{a}}\left({\chi }_{i}^{0}{\chi }_{j}^{0}+{\underline{\chi }}_{i}^{0}{\underline{\chi }}_{j}^{0}\right),$$

(11)

$$=\mathop{\sum}\limits_{a,{\langle ij\rangle }_{a}}2{J}_{a}{\hat{u}}_{{\langle ij\rangle }_{a}}\left({f}_{i}^{{\dagger} }{f}_{j}^{}+{{{\rm{h.c.}}}}\right).$$

(12)

We introduce a hole at site j₀ by switching off all interactions at the corresponding bonds,

$${\tilde{H}}_{d}^{{j}_{0}}={H}_{d}-\mathop{\sum}\limits_{a,j\in {\langle {j}_{0}\,j\rangle }_{a}}2{J}_{a}{\hat{u}}_{{\langle {j}_{0}\,j\rangle }_{a}}\left({f}_{{j}_{0}}^{{\dagger} }{f}_{j}^{}+{{{\rm{h}}.{\rm{c}}.}}\right).$$

(13)

Let $\left\vert {{\Omega }}\right\rangle$ and $\left\vert \tilde{{{\Omega }}}\right\rangle$ 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=| \left\langle \tilde{{{\Omega }}}| {{\Omega }}\right\rangle {| }^{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 the hole Hamiltonian in the eigenbasis at half-filling, ${\tilde{H}}_{d}^{{j}_{0}}={\sum }_{i}{\tilde{\varepsilon }}_{i}\left\vert {\tilde{\alpha }}_{i}\right\rangle \left\langle {\tilde{\alpha }}_{i}\right\vert$,

$$A(\omega )=\mathop{\sum}\limits_{i}\delta (\omega -{\tilde{\varepsilon }}_{i})| \left\langle {\tilde{\alpha }}_{i}| {{\Omega }}\right\rangle {| }^{2}.$$

(14)

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 one-particle excitations on top of the ground state already account for most of the spectral weight⁶⁷, e.g., for J_z/J = 2.5, we find ∫dωA₁(ω) = 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.

Fig. 6: Slow-hole approximation.

The spectral density $D(\omega )=\frac{1}{{L}_{x}{L}_{y}}{\sum }_{{{{\bf{k}}}}}A({{{\bf{k}}}},\omega )$ 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.

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. ¹⁷. 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|\, \ll \,{(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 $\left\vert \tilde{{{\Omega }}}\right\rangle$. 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¹⁷. This can be intuitively understood from the isolated dimer limit J_z/J → ∞, where the ground state consists of dimers on the z-bonds $({\left\vert \uparrow \downarrow \right\rangle }_{z}+{\left\vert \downarrow \uparrow \right\rangle }_{z})/\sqrt{2}$. Inserting either an up- or down-hole on one lattice site will break up the corresponding dimer, e.g., ${\vert \uparrow \circ \rangle }_{z}$. The hole can hop freely back and forth along the same dimer on the z-bond ${\vert \uparrow \circ \rangle }_{z}\to {\vert \circ \uparrow \rangle }_{z}$. Hopping along the x- or y-bonds requires hopping to a neighboring dimer

$$\begin{array}{l}{\vert\uparrow \circ \rangle }_{z}\otimes ({\vert \uparrow \downarrow \rangle }_{\tilde{z}}+{\vert \downarrow \uparrow \rangle }_{\tilde{z}})/\sqrt{2}\\\to ({\vert \uparrow \uparrow \rangle }_{z}\otimes {\vert \circ \downarrow \rangle }_{\tilde{z}}+{\vert \uparrow \downarrow \rangle }_{z}\otimes {\vert \circ \uparrow \rangle }_{\tilde{z}})/\sqrt{2}\end{array}$$

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¹⁷.

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¹⁸, $H(\rho )=\rho {H}_{d}+(1-\rho ){\tilde{H}}_{d}^{{j}_{0}}$; 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,

$${C}_{ij}(\tau )\,\approx \,{e}^{i{E}_{0}\tau }\left\langle {{\Omega }}\right\vert {a}_{i}^{}{{{{\rm{e}}}}}^{-i{H}_{N}\delta }{{{{\rm{e}}}}}^{-i{H}_{N-1}\delta }\ldots {{{{\rm{e}}}}}^{-i{H}_{2}\delta }{{{{\rm{e}}}}}^{-i{H}_{1}\delta }{a}_{j}^{{\dagger} }\left\vert {{\Omega }}\right\rangle ,$$

(15)

where we introduce hole operators ${a}_{j}^{({\dagger} )}$ which locally switch off the couplings. Here, E₀ 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 ${H}_{n}={\sum }_{i}{\varepsilon }_{i}^{(n)}\left\vert {\alpha }_{i}^{(n)}\right\rangle \left\langle {\alpha }_{i}^{(n)}\right\vert$. By insertion of identities after each time step, the correlation function simplifies to a product of overlaps,

$$\begin{array}{rc}{C}_{ij}(\tau )\,\approx \,{e}^{i{E}_{0}\tau }\mathop{\sum}\limits_{{i}_{1},\ldots {i}_{N}}&\left\langle {{\Omega }}\bigg\vert {a}_{i}\bigg\vert {\alpha }_{{i}_{N}}^{(N)}\right\rangle {{{{\rm{e}}}}}^{-i\delta {\varepsilon }_{{i}_{N}}^{(N)}}\left\langle {\alpha }_{{i}_{N}}^{(N)}\bigg\vert {\alpha }_{{i}_{N-1}}^{(N-1)}\right\rangle \\ &\ldots \left\langle {\alpha }_{{i}_{2}}^{(2)}\bigg\vert {\alpha }_{{i}_{1}}^{(1)}\right\rangle {e}^{-i\delta {\varepsilon }_{{i}_{1}}^{(1)}}\left\langle {\alpha }_{{i}_{1}}^{(1)}\bigg\vert {a}_{j}^{{\dagger} }\bigg\vert {{\Omega }}\right\rangle .\end{array}$$

(16)

In practice, again we do not need to sum over all eigenstates $\left\vert {\alpha }_{{i}_{n}}^{(n)}\right\rangle$: 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 $\left\vert {{{\Omega }}}^{(n)}\right\rangle$. We find that indeed the overlaps are very close to one. Thus, we can simplify the time-dependent correlations to

$${C}_{ij}(\tau )\,\approx \,\langle {{\Omega }}\vert {a}_{i}\vert {{{\Omega }}}^{(N)}\rangle \mathop{\prod }\limits_{n=1}^{N-1}({e}^{-i\delta {\varepsilon }_{0}^{(n)}}\langle {{{\Omega }}}^{(n+1)}| {{{\Omega }}}^{(n)}\rangle)\langle {{{\Omega }}}^{(1)}\vert {a}_{j}^{{\dagger} }\vert {{\Omega }}\rangle .$$

(17)

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|\, \ll \,{(J/{J}_{z})}^{4}\,\approx \,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 energy-resolved hole spectrum as measured by ARPES experiments, Eq. (8) or their sum rules. Here, we will consider the spatially averaged local hole spectrum, that is measurable with scanning tunneling microscopy (STM):

$$S(\omega )=\frac{1}{{L}_{x}{L}_{y}}\mathop{\sum}\limits_{\sigma =\uparrow ,\downarrow }\mathop{\sum}\limits_{j}\int{{{\rm{d}}}}\tau {{{{\rm{e}}}}}^{i\tau \omega }\left\langle {\psi }_{0}\right\vert {c}_{j\sigma }^{{\dagger} }(\tau ){c}_{j\sigma }^{}(0)\left\vert {\psi }_{0}\right\rangle .$$

(18)

Although for lattices with a single site per unit cell S(ω) is equal to the spectral density $D(\omega )=\frac{1}{{L}_{x}{L}_{y}}{\sum }_{{{{\bf{k}}}}}A({{{\bf{k}}}},\omega )$, 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.

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).

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 spin-liquid 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³⁹. 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-hole-like 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.

Methods

Physical spins and form factor

In the promising candidate materials exhibiting strong Kitaev interactions, such as 5d iridate compounds Na₂IrO₃ and 4dα-RuCl₃, the partially filled d⁵ orbital is split into e_g and t_2g ($\left\vert xy\right\rangle$, $\left\vert yz\right\rangle$, and $\left\vert zx\right\rangle$) 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 low-energy 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₂IrO₃.

The Kramers doublet, $({c}_{\uparrow }^{{\dagger} },{c}_{\downarrow }^{{\dagger} })$, can be expressed in terms of t_2g orbitals as⁷²

$$\begin{array}{ll}&{c}_{\uparrow }^{{\dagger} }=\frac{1}{\sqrt{3}}\left({d}_{xy,\uparrow }^{{\dagger} }+{d}_{yz,\downarrow }^{{\dagger} }+i{d}_{zx,\downarrow }^{{\dagger} }\right)\\ &{c}_{\downarrow }^{{\dagger} }=\frac{1}{\sqrt{3}}\left(-{d}_{xy,\downarrow }^{{\dagger} }+{d}_{yz,\uparrow }^{{\dagger} }-i{d}_{zx,\uparrow }^{{\dagger} }\right),\end{array}$$

(19)

where ${d}_{m,\sigma }^{{\dagger} }$ 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 ${{{{\bf{d}}}}}_{}^{{\dagger} }=\left({d}_{xy,\uparrow }^{{\dagger} },{d}_{xy,\downarrow }^{{\dagger} },{d}_{yz,\uparrow }^{{\dagger} },{d}_{yz,\downarrow }^{{\dagger} },{d}_{zx,\uparrow }^{{\dagger} },{d}_{zx,\downarrow }^{{\dagger} }\right)$ and define a projector P

$$P=\frac{1}{\sqrt{3}}{\left(\begin{array}{llllll}1&0&0&1&i&\\ 0&-1&1&0&-i&0\end{array}\right)}^{T}.$$

(20)

Then we obtain a compact form of $({c}_{\uparrow }^{{\dagger} },{c}_{\downarrow }^{{\dagger} })={{{{\bf{d}}}}}_{}^{{\dagger} }P$.

The physical spin operators at site j are defined as

$${S}_{j}^{a}=\mathop{\sum}\limits_{m=(xy,yz,zx)}\mathop{\sum}\limits_{s,{s}^{{\prime} }}{d}_{j,m,s}^{{\dagger} }{\sigma }_{s{s}^{{\prime} }}^{a}{d}_{j,m,{s}^{{\prime} }}^{{\dagger} }.$$

(21)

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}^{{\dagger} },...,{c}_{4}^{{\dagger} })$. Using P and P_⊥, one can express the electrons as ${{{{\bf{d}}}}}^{{\dagger} }={{{{\bf{c}}}}}^{{\dagger} }{(P+{P}_{\perp })}^{{\dagger} }$ where ${{{{\bf{c}}}}}^{{\dagger} }\equiv ({c}_{\uparrow }^{{\dagger} },{c}_{\downarrow }^{{\dagger} },{c}_{1}^{{\dagger} },...,{c}_{4}^{{\dagger} })$. 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

$${S}_{j}^{a}\to -\frac{1}{3}\mathop{\sum}\limits_{s,{s}^{{\prime} }}{c}_{s}^{{\dagger} }{\sigma }_{s{s}^{{\prime} }}^{a}{c}_{{s}^{{\prime} }}^{}=-\frac{1}{3}{\sigma }_{j}^{a}.$$

(22)

The tensor of the dynamical hole spectral function is thus defined as

$${{{{\mathcal{A}}}}}_{ms,{m}^{{\prime} }{s}^{{\prime} }}({{{\bf{k}}}},\omega )=\int{{{\rm{d}}}}\tau {e}^{i\tau \omega }\langle {\psi }_{0}| {d}_{{{{\bf{k}}}},ms}^{{\dagger} }(\tau ){d}_{{{{\bf{k}}}},{m}^{{\prime} }{s}^{{\prime} }}^{}(0)| {\psi }_{0}\rangle .$$

(23)

By viewing ${{{\mathcal{A}}}}({{{\bf{k}}}},\omega )$ as a 6 × 6 matrix, it is related to the spectral function for c-fermions by a unitary transformation as ${{{\mathcal{A}}}}({{{\bf{k}}}},\omega )={(P+{P}_{\perp })}^{* }{{{\mathcal{F}}}}({{{\bf{k}}}},\omega ){(P+{P}_{\perp })}^{T}$, where,

$${{{{\mathcal{F}}}}}_{\mu ,\nu }({{{\bf{k}}}},\omega )=\int{{{\rm{d}}}}\tau {e}^{i\tau \omega }\langle {\psi }_{0}| {c}_{{{{\bf{k}}}},\mu }^{{\dagger} }(\tau ){c}_{{{{\bf{k}}}},\nu }^{}(0)| {\psi }_{0}\rangle$$

(24)

with H the Hamiltonian and E₀ 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

$${{{\mathcal{F}}}}\,\approx \left(\begin{array}{ccc}{F}_{\uparrow ,\uparrow }&{F}_{\uparrow ,\downarrow }&0\\ {F}_{\downarrow ,\uparrow }&{F}_{\downarrow ,\downarrow }&0\\ 0&0&{0}_{4\times 4}\end{array}\right).$$

(25)

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}_{\sigma ,{\sigma }^{{\prime} }}$, as

$${{{\mathcal{A}}}}=\frac{1}{3}\left(\begin{array}{cccccc}{F}_{\uparrow ,\uparrow }&-{F}_{\uparrow ,\downarrow }&{F}_{\uparrow ,\downarrow }&{F}_{\uparrow ,\uparrow }&-i{F}_{\uparrow ,\downarrow }&i{F}_{\uparrow ,\uparrow }\\ -{F}_{\downarrow ,\uparrow }&{F}_{\downarrow ,\downarrow }&-{F}_{\downarrow ,\downarrow }&-{F}_{\downarrow ,\uparrow }&i{F}_{\downarrow ,\downarrow }&-i{F}_{\downarrow ,\uparrow }\\ {F}_{\downarrow ,\uparrow }&-{F}_{\downarrow ,\downarrow }&{F}_{\downarrow ,\downarrow }&{F}_{\downarrow ,\uparrow }&-i{F}_{\downarrow ,\downarrow }&i{F}_{\downarrow ,\uparrow }\\ {F}_{\uparrow ,\uparrow }&-{F}_{\uparrow ,\downarrow }&{F}_{\uparrow ,\downarrow }&{F}_{\uparrow ,\uparrow }&-i{F}_{\uparrow ,\downarrow }&i{F}_{\uparrow ,\uparrow }\\ i{F}_{\downarrow ,\uparrow }&-i{F}_{\downarrow ,\downarrow }&i{F}_{\downarrow ,\downarrow }&i{F}_{\downarrow ,\uparrow }&{F}_{\downarrow ,\downarrow }&-{F}_{\downarrow ,\uparrow }\\ -i{F}_{\uparrow ,\uparrow }&i{F}_{\uparrow ,\downarrow }&-i{F}_{\uparrow ,\downarrow }&-i{F}_{\uparrow ,\uparrow }&-{F}_{\uparrow ,\downarrow }&{F}_{\uparrow ,\uparrow }\\ \end{array}\right).$$

(26)

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

$$A({{{\bf{k}}}},\omega )={F}_{\uparrow ,\uparrow }({{{\bf{k}}}},\omega )+{F}_{\downarrow ,\downarrow }({{{\bf{k}}}},\omega ).$$

(27)

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⁷³. Numerical costs grow 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 ${{{\rm{U}}}}(1)\times {{\mathbb{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 x-direction. We favor the flux-free sector explicitly as a ground state in our simulation by adding plaquette terms ∝ − ∑_pW_p. For a periodic or infinite system, it is known that the ground state has to be flux-free⁷⁴.

The time evolution was realized via an MPO representation of the time evolution operator, the W_II operator from ref. ⁷⁵. 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 $\left\vert {\psi }_{0}\right\rangle$ of the pure Kitaev model

$$\left\vert \psi (\tau )\right\rangle ={{{{\rm{e}}}}}^{-iH\tau }{c}_{j\sigma }^{}\left\vert {\psi }_{0}\right\rangle ,$$

(28)

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 not reach the boundaries on the simulated time scales to limit the boundary effects.

To access the spectral function, we employ established MPS methods⁷⁶. Initially, we obtain an MPS approximation for the ground state without hole $\left\vert {\psi }_{0}\right\rangle$ 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

$${C}_{ij}(\tau )=\mathop{\sum}\limits_{\sigma }\left\langle {\psi }_{0}| {{{{\rm{e}}}}}^{i\tau H}{c}_{j\sigma }^{{\dagger} }{{{{\rm{e}}}}}^{-i\tau H}{c}_{i\sigma }^{}| {\psi }_{0}\right\rangle ,\quad i,j\in \{A,B\}.$$

(29)

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 sublattices 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

$${c}_{{{{\bf{k}}}}\sigma }^{}=\frac{1}{\sqrt{N}}\mathop{\sum}\limits_{i\in \{A,B\}}{{{{\rm{e}}}}}^{i{{{\bf{k}}}}{{{{\bf{r}}}}}_{i}}{c}_{i\sigma }^{}.$$

(30)

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

$${c}_{{{{\bf{k}}}},A(B)\sigma }^{}=\frac{1}{\sqrt{N}}\mathop{\sum}\limits_{i\in A(B)}{{{{\rm{e}}}}}^{i{{{\bf{k}}}}{{{{\bf{r}}}}}_{i}}{c}_{i,A(B)\sigma }^{}.$$

(31)

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:

$${c}_{{{{\bf{k}}}}\sigma }^{}=\left\{\begin{array}{ll}{c}_{{{{\bf{k}}}},A\sigma }^{}+{c}_{{{{\bf{k}}}},B\sigma }^{}\quad &{{{\rm{for\,}}}}{{{\bf{k}}}}\in {{{\rm{1.BZ}}}}\\ {c}_{{{{\bf{k}}}},A\sigma }^{}+{{{{\rm{e}}}}}^{i{{{\bf{g}}}}{{{{\bf{r}}}}}_{AB}}{c}_{{{{\bf{k}}}},B\sigma }^{}\quad &{{{\rm{for\,}}}}{{{\bf{k}}}}-{{{\bf{g}}}}\in {{{\rm{1.BZ}}}}\end{array}\right.,$$

(32)

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⁵⁵. In this sublattice spectral function, the lower branch of the free-hole 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⁷⁸ 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.

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.

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.

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.

For the ferromagnetic case, Fig. 9a, 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. 9c, 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. 9d 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 couplings J_z/J = 1 in Fig. 10. This corresponds to the gapless phase of the Kitaev spin liquid⁸. 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.

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.

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. ¹⁹ and consider the splitting of a hole into a charge degree of freedom, the holon${b}_{}^{{\dagger} }$, and a spin degree of freedom, the spinon${f}_{\sigma }^{{\dagger} }$,

$${c}_{i\uparrow }^{}=\frac{1}{\sqrt{2}}\left({b}_{i1}^{{\dagger} }{f}_{i\uparrow }^{}-{b}_{i2}^{{\dagger} }{f}_{i\downarrow }^{{\dagger} }\right),\quad {c}_{i\downarrow }^{}=\frac{1}{\sqrt{2}}\left({b}_{i1}^{{\dagger} }{f}_{i\downarrow }^{}+{b}_{i2}^{{\dagger} }{f}_{i\uparrow }^{{\dagger} }\right).$$

(33)

Here, the holons fulfill the bosonic commutation relations $[{b}_{in}^{},{b}_{jm}^{{\dagger} }]={\delta }_{ij}{\delta }_{nm}$ and the spinons obey fermionic anticommutation relations $\{{f}_{i\sigma }^{},{f}_{j\tilde{\sigma }}^{{\dagger} }\}={\delta }_{ij}{\delta }_{\sigma \tilde{\sigma }}$. By introducing two holon species ${b}_{i1}^{{\dagger} }$, ${b}_{i2}^{{\dagger} }$ we retain the SU(2) gauge redundancy. At each site we must impose the constraint ${f}_{i\uparrow }^{{\dagger} }{f}_{i\uparrow }^{}+{f}_{i\downarrow }^{{\dagger} }{f}_{i\downarrow }^{}+{b}_{i1}^{{\dagger} }{b}_{i1}^{}-{b}_{i2}^{{\dagger} }{b}_{i2}^{}=1$. 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,

$${f}_{i\uparrow }^{}=\frac{1}{\sqrt{2}}\left({\chi }_{i}^{0}+i{\chi }_{i}^{z}\right),\quad {f}_{i\downarrow }^{}=\frac{1}{\sqrt{2}}\left(i{\chi }_{i}^{x}-{\chi }_{i}^{y}\right).$$

(34)

Therefore, we can rewrite the Kitaev model in terms of these spinons:

$${H}_{K}=-\mathop{\sum}\limits_{\alpha ,\sigma }\mathop{\sum}\limits_{{\langle i,j\rangle }_{\alpha }}\left({\tilde{J}}_{\alpha }^{\sigma }{f}_{i,\sigma }^{{\dagger} }{f}_{j,\sigma }^{}+{{{\Delta }}}_{\alpha }^{\sigma }{f}_{i,\sigma }^{{\dagger} }{f}_{j,\sigma }^{{\dagger} }+{{{\rm{h.c.}}}}\right),$$

(35)

where the parameters $\tilde{J}$, Δ are determined self-consistently⁷⁹.

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 decoupling¹⁹

$$\begin{array}{l}{H}_{t}=-t\mathop{\sum}\limits_{\langle i,j\rangle ,\sigma }\left({c}_{i\sigma }^{{\dagger} }{c}_{j\sigma }^{}+{{{\rm{h.c.}}}}\right)\\ \quad\,\,=-\frac{t}{2}\mathop{\sum}\limits_{\langle i,j\rangle ,\sigma }\left({b}_{i1}^{}{b}_{j1}^{{\dagger} }{f}_{i\sigma }^{{\dagger} }{f}_{j\sigma }^{}-\sigma {b}_{i1}^{}{b}_{j2}^{{\dagger} }{f}_{i\sigma }^{{\dagger} }{f}_{j-\sigma }^{{\dagger} }\right.\\ \left.\qquad\,\,-\,\sigma {b}_{i2}^{}{b}_{j1}^{{\dagger} }{f}_{i-\sigma }^{}{f}_{j\sigma }^{}+{b}_{i2}^{}{b}_{j2}^{{\dagger} }{f}_{i-\sigma }^{}{f}_{j-\sigma }^{{\dagger} }+{{{\rm{h.c.}}}}\right)\\ \quad\,\approx -\frac{t}{2}\mathop{\sum}\limits_{\langle i,j\rangle ,\sigma }\left(\langle {b}_{i1}^{}{b}_{j1}^{{\dagger} }\rangle {f}_{i\sigma }^{{\dagger} }{f}_{j\sigma }^{}+{b}_{i1}^{}{b}_{j1}^{{\dagger} }\langle {f}_{i\sigma }^{{\dagger} }{f}_{j\sigma }^{}\rangle \right.\\ \qquad-\,\sigma \langle {b}_{i1}^{}{b}_{j2}^{{\dagger} }\rangle {f}_{i\sigma }^{{\dagger} }{f}_{j-\sigma }^{{\dagger} }-\sigma {b}_{i1}^{}{b}_{j2}^{{\dagger} }\langle {f}_{i\sigma }^{{\dagger} }{f}_{j-\sigma }^{{\dagger} }\rangle \\ \qquad-\,\sigma \langle {b}_{i2}^{}{b}_{j1}^{{\dagger} }\rangle {f}_{i-\sigma }^{}{f}_{j\sigma }^{}-\sigma {b}_{i2}^{}{b}_{j1}^{{\dagger} }\langle {f}_{i-\sigma }^{}{f}_{j\sigma }^{}\rangle \\ \qquad\left.+\,\langle {b}_{i2}^{}{b}_{j2}^{{\dagger} }\rangle {f}_{i-\sigma }^{}{f}_{j-\sigma }^{{\dagger} }+{b}_{i2}^{}{b}_{j2}^{{\dagger} }\langle {f}_{i-\sigma }^{}{f}_{j-\sigma }^{{\dagger} }\rangle +{{{\rm{h.c.}}}}\right)\\ \quad=-\,\frac{t}{2}\mathop{\sum}\limits_{\langle i,j\rangle }\left({W}_{ij}{b}_{i1}^{}{b}_{j1}^{{\dagger} }-{W}_{ij}^{* }{b}_{i2}^{}{b}_{j2}^{{\dagger} }+{{{\rm{h.c.}}}}\right),\end{array}$$

(36)

where we defined ${W}_{ij}={\sum }_{\sigma }\langle {f}_{i\sigma }^{{\dagger} }{f}_{j\sigma }^{}\rangle$ and take the expectation value of the undoped Kitaev spin-liquid ground state. Moreover, we used that $\langle {b}_{i\mu }^{}{b}_{j\mu }^{{\dagger} }\rangle =0$ because after the creation of a holon at site j in the ground state, there cannot be another hole at site i for i ≠ j. Furthermore, we have $\langle {f}_{i\sigma }^{{\dagger} }{f}_{j-\sigma }^{{\dagger} }\rangle =0$ since these pairing terms do not conserve spin parity. Thus, the full Hamiltonian decomposes into H = H_t + H_K, where the first part only acts on the holon and the second part describes the spinon dynamics. Accordingly, for the spectral function, we can employ the spectral building principle similar to the one-dimensional t-J model⁶³. Concretely, the time-dependent correlation function Eq. (29), factorizes into a product of the holon and the spinon part.

$$\begin{array}{rc}{C}_{ij}(\tau )=&\mathop{\sum}\limits_{\sigma }\left(\left\langle {\psi }_{0}\bigg\vert {{{{\rm{e}}}}}^{i{H}_{K}\tau }{f}_{i\sigma }^{{\dagger} }{{{{\rm{e}}}}}^{-i{H}_{K}\tau }{f}_{j\sigma }^{}\bigg\vert {\psi }_{0}\right\rangle \left\langle 0\bigg\vert {{{{\rm{e}}}}}^{i{H}_{t}\tau }{b}_{i1}^{}{{{{\rm{e}}}}}^{-i{H}_{t}\tau }{b}_{j1}^{{\dagger} }\bigg\vert 0\right\rangle \right.\\ &\left.+\left\langle {\psi }_{0}\bigg\vert {{{{\rm{e}}}}}^{i{H}_{K}\tau }{f}_{i\sigma }^{}{{{{\rm{e}}}}}^{-i{H}_{K}\tau }{f}_{j\sigma }^{{\dagger} }\bigg\vert {\psi }_{0}\right\rangle \left\langle 0\bigg\vert {{{{\rm{e}}}}}^{i{H}_{t}\tau }{b}_{i2}^{}{{{{\rm{e}}}}}^{-i{H}_{t}\tau }{b}_{j2}^{{\dagger} }\bigg\vert 0\right\rangle \right)\end{array}$$

(37)

After Fourier transformation to momentum space, the spectrum is given by the convolution of the individual holon and spinon spectral functions A_h(k, ω) and A_sp(k, ω), respectively.

$$\begin{array}{ll}A({{{\bf{k}}}},\omega )\,=\,\int{{{\rm{d}}}}\nu {{{\rm{d}}}}{{{\bf{q}}}}\,{A}_{{{{\rm{sp}}}}}({{{\bf{k}}}}-{{{\bf{q}}}},\omega -\nu ){A}_{{{{\rm{h}}}}}({{{\bf{q}}}},\nu )\\ \qquad\quad\,\,\,\,=\,\mathop{\sum}\limits_{\alpha ,\beta }\int{{{\rm{d}}}}{{{\bf{q}}}}\,{Z}_{{{{\rm{sp}}}}}^{\alpha }({{{\bf{k}}}}-{{{\bf{q}}}}){Z}_{{{{\rm{h}}}}}^{\beta }({{{\bf{q}}}})\delta (\omega -{\varepsilon }_{{{{\rm{sp}}}}}^{\alpha }({{{\bf{k}}}}-{{{\bf{q}}}})-{\varepsilon }_{{{{\rm{h}}}}}^{\beta }({{{\bf{q}}}})),\end{array}$$

(38)

with labels α, β for the multiple bands of the holon and spinon mean-field Hamiltonians corresponding to different dispersions ε(k) and quasiparticle weights Z(k). As discussed in the previous section, the periodicity in ARPES experiments extends to a larger BZ because of the two-site unit cell of the honeycomb lattice. In Eq. (38), the internal integration over q runs only over the first BZ, while k has to be computed for the larger BZ. The dispersion relations are periodic in the 1.BZ, but the quasiparticle weights must be calculated according to Eq. (32).

The convolution above generally leads to a broad distribution of spectral weight. For large t/∣J∣, the charge degree of freedom is assumed to have a considerable dispersion, whereas the spins have possible low energy modes. Therefore, typically, we expect that spinons and 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_ij∣ ≈ 0.23 on x- and y- bonds and ∣W_ij∣ ≈ 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. 11b 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 x- and 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 ${\varepsilon }_{{{{\rm{h}}}}}({{{{\bf{k}}}}}_{{{{\rm{vH}}}}})+{\varepsilon }_{{{{\rm{sp}}}}}^{x,y,z}$. 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.

Fig. 11: Hole spectral functions from MPS and parton mean-field theory.

Spectral function A(k, ω) at low energies for antiferromagnetic J < 0, J_z/J = 2.5 and t/∣J∣ = 3.0 (a) from MPS time evolution and (b) from parton convolution; see Eq. (38). Dashed/ dotted lines correspond to energies of spinons for bond excitations along x- and y-bonds/ z-bonds, respectively, added to the van Hove points of the holons, which are expected to dominate the spectral response.

The MPS spectrum features similar flat bands as well; see Fig. 11a. 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.