Z2 flux binding to higher-spin impurities in the Kitaev spin liquid

Magnetic impurities in strongly correlated electron systems are valuable tools for probing hidden physical phenomena. A well-known example is the Kondo effect, where the scattering of conducting electrons by magnetic impurities in metals or quantum dots reveals crucial insights into the low-energy physics of both the bulk material and the impurities^1,2,3,4,5. This effect, characterized by the screening of the magnetic impurity by the conduction electrons leading to the formation of a Kondo singlet, has profound implications for understanding many-body interactions and has been extensively studied in both theoretical and experimental contexts.

The study of magnetic impurities extends to a variety of systems beyond conventional metals. In topological insulators, for instance, magnetic impurities can break time-reversal symmetry, leading to the opening of a gap at the Dirac point on the surface states and potentially inducing novel magnetic phases^6,7. These systems provide a rich playground for exploring the interplay between magnetism and topology, with potential applications in spintronics and quantum computing.

In low-dimensional spin systems such as quantum spin liquids (QSLs), the introduction of impurities can reveal even more exotic phenomena^8,9,10,11. QSLs, which are characterized by a lack of conventional magnetic order even at zero temperature due to strong quantum fluctuations, offer a unique environment where impurities can induce localized excitations and modify the emergent gauge fields.

In the context of the Kitaev spin liquid (KSL) model¹² -a paradigmatic example of a two-dimensional QSL with fractionalized excitations and emergent gauge fields—introducing impurities, whether magnetic (spin-S sites) or non-magnetic (vacancies)^{13,14,15,16,17,18,19,20,21,22}, can lead to various novel phenomena. These include localized bound states^15,19,23,24, flux binding effects^{13,14,17,18,22}, and modifications in the system’s topological nature^20,22. Such effects not only provide new insights into impurity physics in QSLs but also enhance our understanding of their overall behavior. Due to the localized nature of the low-energy fractionalized excitations in the above phenomena, characterization of KSL based on the signatures in scanning tunneling microscopy (STM) has been proposed in various theoretical works^{21,23,24,25,26,27,28,29,30,31,32,33}. In particular, the above proposals show that the inelastic tunneling spectroscopy can access the real-space spin-spin correlation function of the Kitaev model, which is very sensitive to defects, open edges, and local flux structures.

So far, two types of local impurities have been relatively well-studied in the Kitaev model. The first type is vacancies. It has been demonstrated that vacancies in the Kitaev model lead to almost zero-energy localized bound states and flux-binding effects^{13,14,15,18,19,34}, which can potentially be probed by thermodynamics^35,36 and STM^21,23,24. The second type is spin-S impurities, which are coupled to KSL at a given site via Kondo coupling. The studies of Kondo impurities^16,17,37 have highlighted several remarkable properties of the Kondo effect in the Kitaev model. In the presence of a spin-1/2 Kondo impurity, the fluxes in the three plaquettes adjacent to the impurity site are no longer individually conserved. However, their product (the flux in the impurity plaquette) and all outer fluxes remain conserved¹⁷. Furthermore, a topological transition occurs from the zero-flux state to a bound-flux state attached to the impurity site as a function of Kondo coupling^16,17.

In this work, we theoretically investigate the behavior of S_imp = 1 and S_imp = 3/2 impurities in KSL by means of numerical exact diagonalization (ED) and density matrix renormalization group (DMRG) methods as well as phenomenological models in the Majorana representation. We will focus on the case with a single magnetic impurity, although the case of multiple impurities can be straightforwardly extended.

We show that the behavior of KSL with a magnetic impurity strongly depends on whether the impurity has a half-integer or integer spin. This dependence, which we demonstrate by considering two cases of magnetic impurities, S_imp = 1 and S_imp = 3/2, echoes the recent findings by Ma³⁸, which shows that the nature of the spin-S Z₂ KSL differs based on whether the spin is integer or half-integer. While Ma’s work introduces 2S-flavor Majorana representation for the pure spin-S KSL and identifies the Z₂ gauge fluxes as conserved quantities, we investigate the behavior of a spin-S impurity embedded within a spin-1/2 Kitaev spin liquid, focusing on the impurity magnetization and the ground-state flux sector as functions of varying impurity coupling strength and spin size. This mixed-spin KSL system, though less explored in the literature, holds significant potential for realization in Kitaev materials containing magnetic impurities.

Similarly to vacancies, the magnetic impurities can bind Z₂-fluxes in the lattice. We show that varying the coupling of spin-S impurity with the surrounding spin-1/2 KSL can drive a phase transition between bound-flux and zero-flux sectors. Furthermore, the point at which this transition occurs depends on the magnitude of S_imp. This is significant because, in the presence of a time-reversal symmetry-breaking magnetic field, each Z₂ flux can bind a Majorana zero mode, thereby realizing an Ising anyon governed by non-Abelian statistics¹². Therefore, demonstrating that magnetic impurities can trap Z₂ fluxes provides a pathway to realizing Ising anyons in these systems.

Model

To describe the Kitaev model with magnetic impurities, we begin with the Hamiltonian:

Here, S^μ = σ^μ/2 represents the μ-component of the spin-1/2 operator, with σ denoting the Pauli matrices. ({tilde{S}}^{mu }) represents the μ-component of the spin-S operator of the impurity. The set of impurity sites is denoted by Λ. J > 0 denotes the ferromagnetic coupling strength between spin-1/2 operators away from the impurity, while g > 0 denotes the ferromagnetic coupling strength between the impurity spins and the spin-1/2 operators of the original Kitaev model. We mainly consider the case that a single impurity is not located on the edges.

Since the interaction on each honeycomb bond remains Kitaev-like, we can define a triple-plaquette flux operator in the vicinity of the magnetic impurity as:

where the product is taken over all bonds forming a 12-site plaquette around the impurity site labeled as “0” (see Fig. 1a). This operator captures the flux configuration within the three plaquettes surrounding the magnetic impurity and allows us to analyze the effect of the impurity on the local flux dynamics. A Z₂ flux operator at a plaquette p in the bulk, where there is no impurity spin, is symbolically expressed as ({W}_{p}={2}^{6}{prod }_{jin p}{S}_{j}^{mu }) in the spin-1/2 basis¹². Both W_p and W_I commute with the Hamiltonian H and with each other, taking ± 1 as their eigenvalues, respectively. Consequently, the total Hilbert space is divided into individual flux sector subspaces: ({mathcal{L}},{ = bigoplus }_{{w}_{{p}_{1}},{w}_{{p}_{2}},cdots ,{w}_{I}},{{mathcal{L}}}_{{w}_{{p}_{1}},{w}_{{p}_{2}},cdots ,{w}_{I}}). We refer to the sector with w_I = − 1 as the bound-flux sector and the sector with w_I = + 1 as the zero-flux sector, when all plaquettes not in the vicinity of the impurity have w_p = 1, as introduced in a vacancy case¹³.

a A spin-S impurity embedded in the Kitaev spin liquid. The inset shows site labels around the impurity site (0) and the plaquette labels used in Eq. (5). b, c Typical expectation values of internal flux operators ({{mathbb{W}}}_{a}) including the impurity site (0) in the bound-flux sector (left) and the zero-flux sector (right) for S_imp = 1 and S_imp = 3/2, respectively. d Local magnetization of the impurity spin embedded in spin-1/2 Kitaev spin liquid, calculated via spin-space exact diagonalization in a 24-site cluster.

Full size image

Internal plaquette operators

Around the impurity, the three adjacent plaquette operators must incorporate the higher-spin operator ({tilde{S}}_{0}^{alpha }), such that they are constructed by the unitary operators of π-rotation, as introduced by Baskaran et al.³⁹. We first quote the basic properties of these unitary operators⁴⁰:

where α ≠ β ≠ γ, and (α, β, γ) ∈ (x, y, z) obey cyclic permutations. The above relations are defined on the same site j, otherwise the operators simply commute on different sites. For a general hexagonal plaquette with arbitrary spin on the corners, one can define the plaquette operator as

Specifically, for the three plaquettes shown in Fig. 1a, we can write down (the tilde is dropped for spin-1/2 sites for clarity)

and then it is straightforward to show that ({{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}={(-1)}^{2{S}_{{rm{imp}}}+2{S}_{5}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{1}}). Therefore, one can see that the mutual commutation relation between the internal plaquette operators depends on the size of the spins shared by the two hexagons, which is consistent with the previous work on the mixed-spin Kitaev model⁴¹. In our impurity model, only the impurity site can have a general spin size and all the rest are spin-1/2. This leads to the important properties of the internal plaquette operators:

where a ≠ b and (a, b) ∈ (I₁, I₂, I₃).

Note that a triple-plaquette operator and internal plaquette operators are directly related. They satisfy the following relation:

For a half-integer spin impurity, the conserved value w_I (= ± 1) can be decomposed into the product of ({w}_{{I}_{1}},{w}_{{I}_{2}},) and ({w}_{{I}_{3}}), each taking ± 1. Thus, in the case of the S_imp = 3/2, the triple-plaquette operator W_I and three internal plaquette operators ({{mathbb{W}}}_{a}) are related by ({W}_{I}=-{{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{3}}). There are four possible configurations for (({w}_{{I}_{1}},{w}_{{I}_{2}},{w}_{{I}_{3}})=(+1,+1,+1),(+1,-1,-1),(-1,+1,-1)) and (− 1, − 1, + 1) that can realize the bound-flux sector w_I = −1 (see the left panel of Fig. 1c). Similarly, there are four degenerate internal flux configurations that can realize the zero-flux sector (see the right panel of Fig. 1c). In contrast, for S_imp = 1, the internal ({{mathbb{W}}}_{a}) always takes on a purely imaginary value with some real coefficient (cin {mathbb{R}}). This makes it independent of w_I (see Fig. 1b). These distinct behaviors for different impurity spin values are clearly observed in the DMRG calculations, which we will discuss more later.

Integer/half-integer dependence of the impurity magnetization

Based on the algebra of the internal plaquette operators (6), one can derive an interesting integer/half-integer effect on the magnetization of the impurity spin. Notice that each internal plaquette operator contains only one of the three spin components of the impurity. This leads to the commutation/anticommutation relations between the internal plaquette operators and impurity spin operators:

From Eq. (4), it follows that if the S_imp is a half-integer, the square of the internal plaquette operators is one (({{mathbb{W}}}_{a}^{2}=+1)) and they all mutually commute ((left[{{mathbb{W}}}_{a},{{mathbb{W}}}_{b}right]=0)). Using these properties, we can demonstrate that:

which implies (langle {tilde{S}}_{0}^{z}rangle =0). The derivation is applicable to the other two components, (langle {tilde{S}}_{0}^{x}rangle =langle {tilde{S}}_{0}^{y}rangle =0), because one can always find an internal plaquette operator that anticommutes with the impurity spin component. The same argument applies to the three neighboring components (({S}_{1}^{x}), ({S}_{9}^{y}), and ({S}_{5}^{z})), as well as to all other spin-1/2 operators on the lattice. In contrast, for an integer-spin impurity, ({w}_{{I}_{a}}) takes pure imaginary values, not quantized to Z₂, leading to the possibility of a non-zero local magnetization (langle {tilde{S}}_{0}^{z}rangle, ne, 0). Similarly, the μ-component (μ = x, y, z) of spin-1/2 neighboring with the impurity spin on μ-bond can have a nonzero spin moment; the other two components should be zero.

This spin-size dependence of the impurity magnetization can be confirmed through numerical exact diagonalization, as illustrated in Fig. 1d. Here we perform ED calculation for the spin Hamiltonian in a 24-site cluster on a cylinder geometry with a single impurity to evaluate the z-component of the magnetic impurity’s local spin moment. For any finite g, the local spin moment is numerically zero for both S_imp = 1/2 and 3/2 cases, which aligns precisely with our analytical findings. Conversely, in the S_imp = 1 case, the local moment takes finite values in the range [− 1, 1] at random due to local dynamics induced by the anticommutation relations of internal plaquette operators. The points plotted in Fig. 1d represent the results obtained in a single run. Averaging ({tilde{S}}_{0}^{z}) across multiple independent runs would yield a mean value approaching zero, highlighting a clear qualitative difference from the half-integer impurity case.

Reentrance effect of the bound-flux sector in zero and finite magnetic fields

Here, we focus on the ground-state flux configuration in the presence of the impurity. Using DMRG within the ITensorsGPU.jl package⁴², we numerically determine the ground-state flux configurations as a function of g/J in a 48-site cylinder with a single spin-S impurity located in the bulk [see Fig. 2a]. Figure 2b shows the transitions between the bound-flux and zero-flux sectors. Starting from the bound-flux sector at g = 0, the system undergoes a flux-sector transition under a weak but nonzero coupling g, leading to the zero-flux sector. This behavior is qualitatively consistent across all S_imp cases, indicating the instability of the bound-flux sector in a “quasivacancy” problem as reported in ref.^18,34. However, the transition point g₁ depends on the spin size of the impurity, shifting to the left as S_imp increases. Additionally, for the S_imp = 3/2 impurity case, the system undergoes a second flux-sector transition at the transition point g₂, resulting in the appearance of the reentrant bound-flux sector. Note again that, in both flux sectors of the S_imp = 3/2 impurity case, a triple-plaquette operator W_I and three internal plaquette operators ({{mathbb{W}}}_{a}) are related as ({W}_{I}=-{{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{3}}). Thus, each flux sector is fourfold degenerate.

a The bound-flux sector of a 48-site cluster used in calculations. Numbers of every plaquette in the bulk shows 〈W_p〉, while the number next to the impurity site shows 〈W_I〉. b Flux-sector transition calculated by DMRG with a single spin-S impurity embedded in spin-1/2 Kitaev spin liquid. c–e Phase diagrams of 〈W_I〉 as a function of the magnetic field strength h/J and the impurity-bulk coupling strength g/J for the cases S_imp = 1/2, 1, and 3/2, respectively. We use the same color scale for all three subfigures: the yellow (gray) region corresponds to the bound-flux (zero-flux) sector. White dot lines denote the edges of each pixel.

Full size image

A few remarks are in order. First, the flux gap exhibits a position dependence of the impurity, possibly due to the edge effect on the flux gap as described by Feng et al.⁴³. Specifically, the flux gap decreases in the cylinder geometry as the impurity approaches one of the edges. This positional dependence might contribute to the energy difference between the two flux sectors. For S_imp = 3/2, this effect is mild and only changes the transition points g₁ and g₂. In contrast, for S_imp = 1, we observe a significant qualitative change in the flux-sector transitions in the strong g region due to the position of the impurity, even on the 48-site cylinder. This change is influenced not only by the specifics of the numerical conditions – such as cluster shape, size, and boundary conditions, – but also by the dynamical properties around the impurity, particularly the lack of quantization of internal flux operators ({{mathbb{W}}}_{a}), as described at Eq. (6). Second, the two-impurity case exhibits behavior qualitatively similar to the single-impurity case. Specifically, every impurity spin with S_imp = 3/2 binds the Z₂ flux in a wide parameter range, while for S_imp = 1, there is a strong position dependence on the flux-sector transitions. Third, our main findings are not specific for the 48-site cylindrical cluster. We have confirmed the same trends in other finite-size cylindrical clusters with different shapes. Numerical evidence supporting these arguments can be found in Supplementary information.

We also examined the stability of the bound-flux sector under a uniform magnetic field numerically using DMRG in the same finite-size cluster. For this analysis, we considered the Hamiltonian ({H}_{{rm{total}}}=H+h{sum }_{j,mu }{S}_{j}^{mu }), with the field applied in the [111] direction in the spin basis. We obtained three phase diagrams of 〈W_I〉 as a function of g and h, shown in Fig. 2c–e, corresponding to the cases of S_imp = 1/2, 1, and 3/2, respectively.

In the cases of S_imp = 1/2 and 1, the bound-flux sector at nonzero g is very fragile in the presence of the external field, and this fragility is independent of the impurity’s position (Supplementary Information). In contrast, the reentrant bound-flux sector for a spin-3/2 impurity exhibits some stability for finite field strengths. This stability also depends on the coupling strength g: as g increases, the bound-flux sector tends to withstand stronger fields. This behavior can be understood by noting that the energy difference between the bound-flux and zero-flux sectors under zero magnetic field increases monotonically with g/J in the reentrant flux sector regime. This increase in energy difference provides the bound-flux sector with greater stability in the presence of magnetic fields.

Triple-plaquette analysis: Majorana representation and effective-coupling model

In this section, we aim to understand better the previous numerical findings, including the spin-size-dependent flux-sector transitions and the reentrance of the bound-flux sector. Similar phenomena have been studied in the site-diluted KSL, where a π-flux can be trapped by a true vacancy or quasivacancy^{13,18,21,24,34}. These studies were based on the exact solution of KSL using the Majorana representation for spin-1/2¹². In these studies, some of the authors have shown that the low-energy modes introduced by quasivacancies are highly localized, allowing for a clear distinction in the energy spectrum between zero- and bound-flux sectors^18,24. Recently, on the other hand, the Majorana parton construction was generalized to study uniform spin-S KSLs³⁸. The essential idea of this construction is recognizing spin-S operator as 2S flavors of spin-1/2s with constraints. This motivates us to re-examine our findings in the mixed-spin KSL using the fermionic approach and the multiple-flavor representation, which we outlined in the Methods. To this end, we focus on a triple-plaquette cluster, as shown in the inset of Fig. 1a, to demonstrate that even the system on the minimal cluster can capture the key findings from the previous section.

In the following, we will discuss two fermionic models in the triple-plaquette cluster: the S_imp = 3/2 impurity model with four-body Majorana interaction term, and the effective-coupling model for a general S_imp. In a general case, the Hamiltonian (1) with a magnetic impurity in Majorana fermion representation reads:

where Λ = 0 denotes the impurity site in the triple-plaquette cluster (see, again, the inset of Fig. 1a). Here, we introduce b^μ and c Majorana operators for the bulk spin-1/2 operator and γ-Majorana operators for the impurity spin as used in ref. ³⁸. The Z₂ gauge field on all the 12 edges is conserved and has eigenvalues ({u}_{jk}^{mu }=pm 1). The triple-plaquette flux, w_I, is thus determined by ({w}_{I}=prod {u}_{jk}^{mu }). The second term in Eq. (10) represents the transformed impurity coupling, which becomes a four-Majorana interaction. Since (i{gamma }_{aj}^{mu }{b}_{k}^{mu }) does not commute with H, this quartic term can not be trivially rewritten as a quadratic form. In our ED calculation, this four-Majorana interaction is treated as it is by preparing 2M Majorana fermion operators composed of (M(in {mathbb{N}})) complex fermion operators in a binary number basis.

We visualize in Fig. 3a the triple-plaquette cluster with S_imp = 3/2 in the Majorana representation. In this case, there are 27 Majorana fermions in the whole 13-spin-site system, including 12 from c-Majoranas of spin-1/2s, 12 from γ-Majoranas of the impurity site, and 3 from b-Majoranas at nearest-neighbor sites of the impurity.

a Majorana representation for the spin-3/2 impurity case in the triple-plaquette cluster, following the parton construction by Ma³⁸. The remaining sites are spin-1/2 degrees of freedom, described using Kitaev’s original parton construction¹². b Energy difference, ΔE, between the bound-flux and zero-flux sectors for the S_imp = 3/2 case. Magenta (orange) curve: ΔE computed with (without) projection operators. Green line: w_I calculated in the spin Hamiltonian (1). c ΔE calculated in the effective-coupling model, Eq (18). The coefficients are given in the main text. From S_imp = 1/2 to higher impurity spins, the substitution (Delta to sqrt{2{S}_{{rm{imp}}}}Delta) is applied.

Full size image

Using ED for the many-body Majorana Hamiltonian (10) under projections for the spin-3/2 impurity, we calculate the energy difference between the bound-flux (w_I = − 1) and zero-flux (w_I = 1) states of the triple-plaquette cluster. This energy difference, defined as

is evaluated as a function of g/J (see the magenta curve in Fig. 3b). When ΔE < 0 (ΔE > 0), the system realizes the bound-flux (zero-flux) sector at a given g/J, which is illustrated as a yellow (gray) shaded area in the figure. We observe two key behaviors: (i) the initial transition from the bound-flux sector to the zero-flux sector at the transition point g₁, and (ii) the reentrant transition back to the bound-flux sector at the second transition point g₂. For g > g₂, ΔE decreases monotonically, indicating the stability of the reentrant bound-flux sector. Moreover, these two transition points match exactly with the flux-sector transition points obtained in ED of the spin-basis Hamiltonian. The calculated w_I in spin-basis ED is shown as the green line in the Fig. 3b. It is also worth mentioning here that the correct behavior of flux-sector transition calculated in the Majorana representation is only obtained under the presence of projection operators P_a and ({P}_{tilde{S} = 3/2}) (see Methods for details). Without these projection operators, as shown by the orange curve in Fig. 3b), the reentrant bound-flux sector does not appear, emphasizing the need for proper constraints on the four-body Majorana terms.

The presented numerical results confirm that even a minimal cluster with a S_imp = 3/2 impurity can qualitatively capture the flux-sector transitions and the reentrant effect, provided a proper treatment of projection. It reinforces the localized picture of the flux-binding effect by a site defect.

To further extract the essential ingredient that leads to a reentrant bound-flux sector, we propose a heuristic approach based on the effective coupling between the c-Majorana eigenmodes of the 12-site plaquette and the Majorana zero modes (({gamma }_{a}^{0})) introduced by the impurity. In the first constituent, the L = 12 plaquette is considered as a fermion hopping problem on a ring, where L denotes the number of sites on the ring⁴⁴. The difference between the zero- and π-flux sectors is translated into the difference between the periodic and antiperiodic boundary conditions (PBC and APBC) of the ring (see Fig. 4a). Here, the π-flux sector corresponds to the bound-flux sector we have discussed. If the nearest-neighbor hopping strength is J/4 (corresponds to the Kitaev coupling in the spin Hamiltonian), the free-fermion Hamiltonian on the ring reads:

with the corresponding energy eigenvalues

where n is the integers from 0 to L − 1. Since the ground-state energy of the system is calculated by the sum of all negative-energy modes, it is straightforward to conclude that L = 4n + 2 favors the zero-flux sector (periodic boundary) and L = 4n favors the π-flux sector (antiperiodic boundary). This is consistent with the prediction by Lieb’s theorem of flux configuration⁴⁵, even only a single plaquette being considered here⁴⁶. Note that in this L = 4n case, ϵ_PBC contains two zero modes while ϵ_APBC contains no zero modes under zero magnetic fields. This is the key difference that helps us understand the flux-sector transition in the presence of impurity.

a The energy spectrum of the fermionic ring with periodic or antiperiodic boundary. The former corresponds to the zero-flux sector of W_I, while the latter corresponds to the bound-flux sector. b The eigenfunctions for some selected eigenmodes. The sites in red can directly hop to the additional impurity fermion. However, if the sum over amplitudes on the red sites vanishes, the effective coupling will be zero.

Full size image

To tackle the quasivacancy or impurity problem with S_imp = 1/2, we add an additional fermion that only directly couples to three sites of the L = 12 ring with effective coupling strength Δ. This Δ can vary for different eigenmodes of the ring coupled to the impurity fermions, but in general, it is proportional to g. Therefore, the effective-coupling Hamiltonian is a tight-binding matrix between some of the eigenmodes of the ring and the impurity fermions. By the symmetry of the wavefunction, only one zero mode (ϵ = 0) and one particle-hole pair (ϵ = ± αJ) of the periodic ring (i.e., zero-flux sector) can hybridize with the impurity (see Methods). This leads to a simple tight-binding matrix:

where Δ_i ~ g is the effective coupling strength and we assume Δ_i ≡ Δ for simplicity. This Hamiltonian gives the eigenvalues of the effective-coupling model when Δ ≪ 1:

where the prime is added to effective-coupling model eigenvalues ({epsilon }^{{prime} }) in order to be distinguished from the hopping-ring model eigenvalues ϵ. One can simply understand the above results by the perturbation theory. If the eigenmode of the L = 12 ring is a zero-energy mode, the effective coupling results in a degenerate perturbation theory with energy correction Δ. On the other hand, if the eigenmodes of the ring have finite energy, the non-degenerate perturbation theory gives rise to second-order corrections. Therefore, we can conclude that the ground-state total energy, which is the sum of all negative eigenvalues, has the general expressions (in the unit of J):

where A₁, B₁, and C₁ are positive constants. In the bound-flux case, there is no zero-energy eigenmode in the plaquette model, so the ground-state energy expression does not include the linear term in Δ when Δ ≪ 1:

In the effective-coupling model, the ground-state energy difference defined in Eq. (11) translates into (Delta E={E}_{{rm{APBC}}}^{{rm{eff}}}-{E}_{{rm{PBC}}}^{{rm{eff}}}), which is the measure of the flux-sector transition. It can be modeled as

and the coefficients are fitted by the exact-diagonalization result for the 13-site fermion-hopping model, which gives A₁ ≈ 1.163, (B₁ − B₂) ≈ − 1.079, and (C₁ − C₂) ≈ − 0.261. Apparently, one can see that when Δ ≪ 1, ΔE is a concave quadratic function and provides the possibility of first (bound-to-zero) and second (zero-to-bound) flux-sector transitions. However, if the predicted second transition happens at the strength ({Delta }_{2}^{* }) beyond the validity of the quadratic approximation, it may not be seen in the exact diagonalization result. This is exactly what happens for the S_imp = 1/2 case shown in the Fig. 3c, where the exact diagonalization curve starts to deviate from the quadratic curve at Δ ≈ 0.5J.

The next important question is how to incorporate the higher-spin impurity. Here, we conjecture that the essential ingredient is the additional ({gamma }_{a}^{0}) Majorana fermions that provide more entries of the tight-binding matrix. For example, for S_imp = 3/2, one introduces three additional zero modes and considers the perturbative effects on the zero mode (ϵ₀ = 0) and the finite-energy modes (ϵ_α = ± αJ) separately:

This indicates that the higher-spin effect can be incorporated by making the substitution (Delta to sqrt{2{S}_{{rm{imp}}}}Delta). In Fig. 3c, we show that the higher-spin impurity in the effective-coupling model simply shifts the quadratic curve to the left. This implies that for the higher-spin impurity case, the first transition point ({Delta }_{1}^{* }) is smaller, and the second transition ({Delta }_{2}^{* }) is more likely to happen when Δ ≪ 1 is still valid. This observation from the effective-coupling model, despite its over-simplicity, is consistent with the DMRG results as shown in Fig. 2a. Based on our findings, one may expect that higher-spin impurity cases such as spin-2, 5/2, and so on, also tend to bind the W_I flux under the (Delta to sqrt{2{S}_{{rm{imp}}}}Delta) substitution. However, it is important to note that our effective model does not account for the effect of ({gamma }_{a}^{mu }), which may lead to unexpected results for even larger spins. Therefore, the precise nature of even higher-spin cases requires further numerical evidence to support the stability of the bound-flux sector.

In this study, we investigated the behavior of S_imp = 1 and S_imp = 3/2 impurities in KSL, focusing on their effects on the flux-sector transitions and the stability of the bound-flux sector. Our analysis, using ED and DMRG methods, along with the phenomenological model in the Majorana representation, revealed several key findings: First, the local behavior of KSL around a magnetic impurity strongly depends on whether the impurity has a half-integer or integer spin. This dependence was demonstrated by considering the cases of S_imp = 1 and S_imp = 3/2, and analyzing their impact on flux-sector transitions. Second, magnetic impurities can bind Z₂ fluxes in the lattice, similar to vacancies and quasivacancies. We observed a phase transition between bound-flux and zero-flux sectors by varying the impurity coupling strength, with the transition point dependent on the impurity’s spin magnitude. Third, for S_imp = 3/2 impurities, a reentrant bound-flux sector was observed, remaining stable under finite magnetic fields. This stability increases with the coupling strength g. Fourth, our ED calculations for fermionic Hamiltonians in the triple-plaquette cluster revealed the energy differences between the bound-flux and zero-flux states, which are in good agreement with the flux-sector transition points identified in the DMRG study of the spin-basis Hamiltonian. Proper constraints by projection operators are essential for accurate flux-sector transition behavior.

The ability of spin-S impurities to bind Z₂ fluxes and its stability against the external magnetic field has practical implications. When time-reversal symmetry is broken, flux binding to the impurity site results in the formation of localized Majorana zero modes, which are essential for realizing Ising anyons. These anyons exhibit non-Abelian statistics, making them valuable for topological quantum computation.

Compared to an Ising anyon bound at a vacancy¹³, the anyon found at the magnetic impurity site offers a more advantageous way to access low-energy Majorana-bound states through its magnetic channel. This unique feature of entangled Ising anyons can be observed in the dynamical correlation function of impurity spins, especially by focusing on the low-energy spectra.

Spectroscopy techniques with high spatial resolution, such as STM or nitrogen-vacancy (NV) center magnetometry, could further advance our understanding of the ground-state flux sector in Kitaev spin liquid phases. The low-energy spectrum carries crucial information about the ground-state flux sector, and both STM and NV center magnetometry are capable of probing these spectral features. STM might reveal the dynamical correlation function of impurity spins through spin-dependent tunneling^21,23,24,30, while NV center magnetometry could detect magnetic noise linked to the same correlation function⁴⁷ as well as emergent gauge field⁴⁸. Observations of low-energy features, including Majorana zero modes, could potentially inform efforts toward realizing topological qubits composed of Ising anyons.

In addition, the even-odd effect on magnetization is another potential signature worth investigating. Local measurements of spin moments, such as those achievable via NV center magnetometry in the absence of a magnetic field, could provide clear evidence of this effect for a magnetic impurity embedded in Kitaev spin liquid phases. However, the spin Hamiltonian discussed in this work may be too simplified to accurately capture magnetic properties in real candidate materials¹¹. Furthermore, introducing more than two magnetic impurities into the system could result in an RKKY-like interaction mediated by itinerant Majorana fermions in the bulk^37,42,49. While this interaction could obscure or destabilize the even-odd effect, it may simultaneously promote magnetic ordering of impurity magnetic moments, even within an otherwise spin-liquid phase.

Implementation of DMRG

All DMRG calculations were performed using the NVIDIA Data Center GPU R470 Driver with the ITensorsGPU.jl package⁴². To ensure qualitative accuracy, determined by the expectation value of all plaquettes for both W_p in the bulk and W_I (plus, even ({{mathbb{W}}}_{a}) for a half-integer impurity case) at impurity sites, we required an adequate bond dimension d depending on the impurity size and a good energy tolerance δE ≤ 1 × 10⁻⁷ while satisfying the cutoff at each sweep to be ≤ 1 × 10⁻⁹. For instance, in the single spin-3/2 impurity problem, we set the maximum bond dimension ({d}_{max }) to 3000 to ensure the cutoff condition. Additionally, we performed DMRG calculations five times independently for each parameter point and selected the result with the lowest ground state energy realizing reasonable flux expectation values.

Construction of projection operators for higher-spin

Here, we first review the Majorana representation for arbitrary spin size introduced by Ma³⁸, and then construct projection operators applied in the ED calculation of the many-body Majorana Hamiltonian for the S_imp = 3/2 case. The starting point is to consider a spin-(tilde{S}) operator as (2tilde{S}) of spin-1/2s:

where S^μ denotes the μ-component of spin-1/2 operator with σ being Pauli matrices, and a indicates the flavor degree of freedom. Then, the Majorana representation is applied for each spin-1/2 as ({sigma }_{a}^{mu }=i{gamma }_{a}^{mu }{gamma }_{a}^{0}) (The order of Majorana operators is slightly modified from the original to follow Kitaev’s Majorana representation¹².). Here four Majorana fermions γ^x, γ^y, γ^z and γ⁰ are introduced for each flavor. As a result,

Note that for spin-1/2 operators, we use the usual symbol of Majorana fermions b^x, b^y, b^z and c instead of ({gamma }_{a}^{eta }) (η = x, y, z, 0). The local Z₂ gauge field operators ({u}_{jk}^{mu }=i{b}_{j}^{mu }{b}_{k}^{mu }), which connect two spin-1/2s at sites j and k on the μ-bond, commute with the Hamiltonian H.

Since the Hilbert space is expanded in this representation- both for the bulk spin-1/2s and impurity sites-we need two kinds of projection operators to ensure the correct physical states are selected. The first one is required to ensure the commutation relations of Pauli matrices at each spin-1/2 or flavor in the Majorana representation, denoted by the condition ({C}_{gamma }^{a}:{D}_{a}={gamma }_{a}^{x}{gamma }_{a}^{y}{gamma }_{a}^{z}{gamma }_{a}^{0}=1) in ref. ³⁸. Resulting operator is

Note that P_a satisfies the condition of a projection operator automatically as ({({P}_{a})}^{2}={P}_{a}).

The second projection operator, which is required only for higher-spin ((tilde{S} ,>, 1/2)) cases, mixes different flavors to ensure (| {tilde{S}}^{2}| =tilde{S}(tilde{S}+1)). This condition is represented as ({C}_{s}:{sum }_{mu }{(mathop{sum }nolimits_{a = 1}^{2tilde{S}}{gamma }_{a}^{mu }{gamma }_{a}^{0})}^{2}=-4tilde{S}(tilde{S}+1)) in ref. ³⁸. For (tilde{S}=3/2), corresponding operator can be represented as

It is worth mentioning that P_a is necessary for ({P}_{tilde{S} = 3/2}) to satisfy the condition of projection operator since

In the Majorana Hamiltonian with four-body interactions in the vicinity of the spin-3/2 impurity, these two kinds of projection operators ({P}_{tilde{S} = 3/2}) and P_a for a = 1, 2, 3 are essential for accurately calculating the energy difference between the two flux sectors.

Eigenenergies and eigenfunctions of the periodic ring

Here we provide details of the L = 12 periodic and antiperiodic rings used in the effective-coupling model. The energy-level spectrum calculated by Eq. (13) is shown in Fig. 4a. This spectrum can be verified by the diagonalization of the tight-binding model as well. The periodic ring, which corresponds to the zero-flux sector, contains two zero-energy modes. This is the crucial difference from the antiperiodic ring, which corresponds to the bound-flux sector. It is important to note that not all eigenmodes can couple to the impurity. First, the impurity fermion only directly hops to three nearest-neighbor sites, which are colored in red in Fig. 4b. Second, one can straightforwardly symmetrize the eigenmodes for each degenerate pair of modes, such that only some of the modes couple to the impurity fermion. Specifically, the 12-site model on a periodic ring has a C_3v point-group symmetry, such that only three modes are subject to the transformation of the trivial irreducible representation and couple to the impurity site. Because the hopping integral is based on the overlap of one site (impurity) and the three neighboring sites simultaneously, one can simply assume ψ_imp = 1 as the impurity wavefunction and the sum of the amplitude of the three sites determines the effective coupling. Specifically, for a given eigenmode of the fermionic ring ψ_ring, if

with j ∈ (3, 7, 11) in Fig. 4b, the eigenmode is simply decoupled to the impurity site. From Fig. 4b, we see that only the eigenmodes with ϵ = ± J/2 and one of the zero modes ϵ = 0 can couple to the impurity. This is numerically verified by the exact diagonalization of the 13-site tight-binding matrix. A similar analysis can be made on the antiperiodic ring for the bound-flux sector. However, since there is no zero-energy mode in the spectrum, the perturbative effect on the eigenmodes by the impurity fermion starts from the second-order correction. This leads to the qualitative argument based on our effective-coupling model, which is discussed in the main text.