Electric polarization near vortices in the extended Kitaev model

We formulate a Majorana mean-field theory for the extended $JK\Gamma$ Kitaev model in a magnetic Zeeman field of arbitrary direction, and apply it for studying spatially inhomogeneous states harboring vortices. This mean-field theory is exact in the pure Kitaev limit and captures the essential physics throughout the Kitaev spin liquid phase. We determine the charge profile around vortices and the corresponding quadrupole tensor. The quadrupole-quadrupole interaction between distant vortices is shown to be either repulsive or attractive, depending on parameters. We predict that electrically biased scanning probe tips enable the creation of vortices at preselected positions. Our results open new perspectives for the electric manipulation of Ising anyons in Kitaev spin liquids.


INTRODUCTION
A hallmark of Kitaev spin liquids is the fractionalization of spin-1/2 local moments into Majorana fermions and a Z 2 gauge field [1][2][3][4][5][6][7][8][9].When time reversal symmetry is broken by an external magnetic field, both types of excitations become gapped, and vortices of the Z 2 gauge field bind Majorana zero modes that behave as non-Abelian anyons.These properties can be demonstrated in the exactly solvable Kitaev honeycomb model [1].Since the observation that the bond-directional exchange interactions of the pure Kitaev model are realized in quasi-two-dimensional Mott insulators with strong spinorbit coupling [10], identifying signatures of fractional excitations in Kitaev materials has become a major goal of condensed matter physics [11][12][13][14].Most notably, there is evidence for a half-quantized thermal Hall conductance in the candidate material α-RuCl 3 at intermediate temperatures and magnetic fields, but its interpretation in terms of chiral Majorana edge modes remains controversial [15][16][17][18].This ambiguity calls for alternative experimental probes that may help distinguish a Kitaev spin liquid from a more conventional partially polarized phase with topological magnons [19,20].
A promising route to detect and manipulate the fractional excitations of Kitaev spin liquids is to exploit their nontrivial responses to electrical probes.Theoretical proposals in this direction include electric dipole contributions to the subgap optical conductivity [21,22], scanning tunneling spectroscopy [23][24][25][26][27], interferometry in electrical conductance [28,29], and electric polarization and orbital currents associated with localized excitations [30,31].In fact, the charge polarization in Mott insulators can be captured by an effective density operator written in terms of spin correlations in the low-energy sector [32,33].The effective density operator for Kitaev materials was derived in Ref. [30] starting from the multiorbital Hubbard-Kanamori model in the ideal limit where the dominant exchange path only generates the pure Kitaev interaction [10].The electric field effects then work both ways.On the one hand, the inhomogeneous spin correlations around a Z 2 vortex imply that vortices pro-duce an intrinsic electric charge distribution.On the other hand, vortices are attracted by electrostatic potentials that locally modify exchange couplings, and this effect can be used to trap and move anyons adiabatically [30,34].
In this work we generalize the theory of the electric charge response in Ref. [30] to consider the generic spin model for Kitaev materials [35,36].Our starting point is the three-orbital Hubbard-Kanamori model which takes into account sub-dominant hopping processes that, in addition to Kitaev (K) interactions, also generate Heisenberg (J) and off-diagonal (Γ) exchange interactions.Using perturbation theory to leading order in the hopping parameters, we derive an expression for the effective charge density operator in the Mott insulating phase that contains all two-spin terms allowed by symmetry.Since the additional interactions spoil the integrability of the pure Kitaev model, we compute spin correlations using a Majorana mean-field theory.This type of approximation has been applied to map out the ground state phase diagram and to compute response functions of the extended Kitaev model [37][38][39][40][41][42][43][44][45][46].Here we generalize the mean-field approach to treat position-dependent order parameters in the case where translation symmetry is broken by the presence of vortices in the Z 2 flux configuration.Including a Zeeman coupling, we show that the spatial anisotropy of the charge distribution around a vortex varies with the direction of the magnetic field and can be quantified by the components of the electric quadrupole moment.We also discuss how a local electrostatic potential renormalizes the couplings in the extended Kitaev model and gives rise to an effective attractive potential for vortices.Remarkably, the effect is stronger in the presence of non-Kitaev interactions, and we find that it is possible to close the vortex gap by means of electric modulation of the local spin interactions.
The local degrees of freedom of Kitaev materials are transition metal ions with 4d 5 or 5d 5 electronic configuration and strong spin-orbit coupling [4,5].In the presence of the crystal field of an octahedral ligand cage, this configuration is equivalent to a single hole in a t 2g orbital.Starting from a three-orbital Hubbard-Kanamori Hamiltonian on the honeycomb lattice, in the presence of a Zeeman coupling to an external magnetic field h, we find from a projection scheme that the low-energy effective spin Hamiltonian is given by the extended Kitaev (aka JKΓ) model [35], where σ i denotes the vector of the pseudospin-1/2 Pauli operators at site i.Moreover, i and j are nearest neighbors, J ij is the bond-dependent exchange matrix, and the indices α, β, γ ∈ {x, y, z} = {1, 2, 3} label both spin components and bonds on the honeycomb lattice.We denote by ⟨ij⟩ γ a nearest-neighbor bond of type γ with site i belonging to sublattice A and j to sublattice B. For bond ⟨ij⟩ z , we have The exchange matrices for x and y bonds follow by cyclic permutation of the spin and bond indices.The ideal Kitaev case with J = Γ = 0 corresponds to a single hopping path mediated by ligands on edge-sharing octahedra with ideal 90 • bonds [10].Numerical studies show that the Kitaev spin liquid phase is stable in the regime |Γ|, |J| ≪ |K| [35,[47][48][49].For estimates of the hopping and exchange parameters for α-RuCl 3 , see for instance Refs.[5,50].In this material, one finds a ferromagnetic Kitaev coupling (K < 0) and the leading perturbation to the idealized Kitaev model is given by 0 < Γ < |K|.
We employ a mean-field approximation for calculating spin correlations in the extended Kitaev model and to verify the stability of the spin liquid phase against integrability-breaking perturbations.For J = Γ = h = 0, the model can be solved exactly [1] using the Kitaev representation σ γ i = ic 0 i c γ i in terms of four Majorana fermions which obey (c µ i ) † = c µ i and {c µ i , c ν j } = 2δ ij δ µν .Throughout, we use indices µ, ν, ρ ∈ {0, 1, 2, 3} to denote all four fermion flavors, in contrast with α, β, γ ∈ {1, 2, 3}.Physical states must respect the local constraint The algebra of the spin operators can be satisfied using different representations [51].It is convenient to write the Kitaev representation in terms of the vector c i = (c 0 i , c 1 i , c 2 i , c 3 i ) T and the antisymmetric matrices N γ defined by Instead of imposing D i = +1, we use the equivalent constraint [52] c Note that the constraints c T i G γ c i = 0 for γ = x, y, z are redundant.If the constraint is implemented exactly, it suffices to impose it for a single value of γ.However, when treating the constraints (3) numerically through the corresponding Lagrange multipliers λ γ i [42,44], it is advantageous to enforce them in a symmetric manner for all three values of γ.We thereby rewrite the spin Hamiltonian as We decouple the quartic terms using two types of realvalued mean-field parameters, which obey For the exactly solvable Kitaev model, one finds that U µν ij is diagonal in the indices µ, ν.In particular, the components U γγ ij are related to the static Z 2 gauge field and take values U γγ ij = ±1 when i, j form a nearest-neighbor γ bond, and U γγ ij = 0 otherwise.Thus, U γγ ij can be viewed as an "order parameter" for the Kitaev spin liquid phase.For comparison with the exact solution, we also define W p = ⟨ij⟩γ ∈p U γγ ij , where p is a hexagonal plaquette.In the pure Kitaev model, W p is identified with the gaugeinvariant Z 2 flux, and the ground state lies in the sector with W p = +1 for all plaquettes.States with W p = −1 at isolated plaquettes are associated with vortex excitations [1].Besides the link variables U µν ij , in the meanfield approach we also consider the on-site fermion bilinears V µν i .It follows from the Kitaev representation that V 0γ i = ⟨σ γ i ⟩.Moreover, the constraint in Eq. ( 3) implies Thus, there are only three independent components of V µν i at each site, and they are related to the local magnetization induced by the external magnetic field.In the limit |h| ≫ |K|, |J|, |Γ|, we expect to encounter a partially polarized phase characterized by V µν i ̸ = 0 while U µν ij = 0 for all bonds.For further detail, see the Methods section.

Homogeneous case
We first describe the mean-field solution for the homogeneous case, i.e., in the absence of vortices.If the ground state does not break spin rotation or lattice symmetries, as in the Kitaev spin liquid phase, the matrices U µν ij depend only on the bond type γ, and we set U µν ij = U µν γ for bonds ⟨ij⟩ γ .Moreover, V µν i = V µν becomes a constant matrix.More generally, we can allow these parameters to vary with the sublattice within larger unit cells to describe magnetically ordered phases.We then solve the mean field self-consistency equations using a Fourier transform of the Majorana modes in the thermodynamic limit.As a first step, we have verified that our mean-field approach recovers the exact results for the Kitaev model [1] when we set Γ = J = h = 0.The resulting dispersion relation of Majorana fermions is depicted by dashed lines in Fig. 1.In this case, the only dispersive band is associated with the fermion c 0 .This band is gapless with a Dirac spectrum near the K point in the Brillouin zone (BZ).In addition, there are three degenerate flat bands associated with the fermions c γ , which are related to the static gauge variables U γγ γ (whose value is independent of γ).
Moving away from the exactly solvable point, we find that all bands become dispersive.For h = 0 and K, J, Γ ̸ = 0, our results are in quantitative agreement with a previous mean-field calculation [37].Our approach also allows us to take into account the magnetic field nonperturbatively.Figure 1 shows the dispersion for a magnetic field pointing along the crystallographic c direction (perpendicular to the honeycomb plane), with unit vector ĉ = 1 √ 3 (1, 1, 1).Here the coordinates are specified in terms of the crystallographic axes x, ŷ and ẑ of the ligand octahedra.For later reference, the in-plane unit vectors are â = 1 √ 6 (1, 1, −2) and b = 1 √ 2 (−1, 1, 0).As shown in Fig. 2, the magnetic field opens up a gap in the fermion spectrum, as expected for the non-Abelian Kitaev spin liquid phase.As we increase the magnetic field, the gap at the K point increases, but the gap at the Γ point decreases.The fermion gap ∆ f is given by the minimum between the energies at the K and Γ points in the BZ.If these energies cross, ∆ f exhibits a kink at the corresponding value of h (e.g., for Γ = 0 in Fig. 2).As we increase the magnetic field, we encounter a critical value h c at which the gap either changes discontinuously, as in a first-order transition (e.g., for Γ = −0.1|K| in Fig. 2), or it vanishes and varies continuously across the phase transition (e.g., for Γ = 0.1|K| in Fig. 2).For h = hĉ and h ≪ h c , the fermion gap increases with the magnetic field as ∆ f ∝ h 3 , as expected from perturbation theory [1]; see the inset in Fig. 2. For general field directions, the fermion gap behaves as ∆ f ∝ h x h y h z , closing when one component of h vanishes.
We further assess the stability of the Kitaev spin liquid phase by evaluating the Z 2 flux parameter.In a homogeneous ground state, we have W p = (U γγ γ ) 6 .The result for the extended Kitaev model with J = 0 and Γ, h ̸ = 0 is shown in Fig. 3 for a magnetic field along the ĉ direction and for an in-plane field along the â direction (perpendicular to the z bonds).As expected, U γγ γ decreases as we increase h or Γ.The dots in this figure mark the transition where the gap ∆ f vanishes continuously.Note that U γγ γ varies smoothly across the continuous transition for h ∥ ĉ and Γ > 0.
The results in Figs. 2 and 3 allow us to determine the parameter regime where both U γγ γ and ∆ f vary smoothly and take values comparable to those at the exactly solvable point.In this regime, we expect the mean-field approach to yield qualitatively correct results for the charge response of the Kitaev spin liquid phase.By contrast, the regime of strong magnetic fields should be identified with the partially polarized phase, whereas the regime of large |Γ| or |J| harbors magnetically ordered phases [ 35,44,48,53].Here we do not explore the various phases of the extended Kitaev model, whose nature is not completely settled [36].Nevertheless, our mean-field results reproduce qualitative features of phase diagrams reported in the literature.For instance, we find that adding Γ > 0 increases the critical magnetic field along the ĉ direction, but the Kitaev spin liquid phase shrinks as we tilt the field towards the plane, in agreement with exact diagonalization results [48].However, in general the mean-field approach overestimates the value of the critical magnetic field for a ferromagnetic Kitaev coupling in comparison with more accurate numerical meth- The ratio between the exchange couplings calculated using Eq. ( 22) are Γ/|K| = 0.20 and J/|K| = −0.02.We set the magnetic field h/|K| = 0.2ĉ.The solid line marks the zigzag path considered in Fig. 5.

Vortex charge density profile
Inhomogeneous spin correlations can bring on a charge redistribution in Mott insulators [32,33].We here discuss the charge density profile induced by the presence of Z 2 vortices in a Kitaev spin liquid.In the Methods section, we derive the effective charge imbalance operator in terms of two-spin operators and show how to compute its expectation value ⟨δn l ⟩ at lattice site l using the Majorana mean-field approach.
We consider an inhomogeneous state in which translation symmetry is broken by the presence of vortices.
In this case, we analyze the mean-field Hamiltonian on a finite system with linear size L along the directions of the primitive lattice vectors ê1 = 2 b, imposing periodic boundary conditions.To create vortices, we initialize the mean-field parameters in a configuration where we flip the sign of U µν ij on bonds crossed by open strings.In the pure Kitaev model, this procedure generates exact eigenstates with two localized vortices at the ends of the string.In the extended Kitaev model, vortices become mobile excitations with effective bandwidths governed by the integrabilitybreaking perturbations [57,58].In fact, for sufficiently large values of these perturbations, near the border of the Kitaev spin liquid phase in Fig. 3, we observe that the vortex positions vary as we iterate the self-consistency equations.When this happens, the string length decreases and the vortices move closer to each other un- til they annihilate, and the mean-field solution converges to the vortex-free ground-state configuration.However, for |Γ|, |J|, |h| ≪ |K| and well separated vortices, we find a self-consistent solution with (metastable) localized vortices which corresponds to a local energy minimum in this sector of the Hilbert space.These results seem consistent with the real-time dynamics described by time-dependent mean-field theory, which show that only when the perturbations are strong enough do vortices become mobile as signaled by the time decay of the fermion Green's function [46].In reality, the lifetime of a vortex is limited by processes in which two vortices meet and annihilate [58], and can become arbitrarily long at low temperatures due to the low vortex density, see the Supplemental Material (SM) [59].Focusing on the regime of small perturbations, we can then compute static spin correlations near vortices using position-dependent mean-field parameters U µν ij and V µν i .We consider a configuration with four equally spaced vortices, see inset of Fig. 5(b), which preserves rotational symmetries and minimizes finite-size effects as compared to a two-vortex configuration.Unless stated otherwise, we use L = 40, so the distance between vortices is 20 unit cells.The charge imbalance near a vortex is then effectively the property of a single vortex and finite-size effects only appear in long-distance tails [59].
In Ref. [30], the charge imbalance profile in the vicinity of a vortex was investigated within the exactly solvable

Kitaev model [1]
setting K γ = K for isotropic Kitaev interactions.The three-spin interaction breaks time-reversal symmetry while preserving integrability.The coupling constant derived from perturbation theory in the magnetic field is [1] where ∆ 2v ≈ 0.263|K| [60] is the energy gap for creating two adjacent vortices at zero magnetic field.The prefactor in Eq. ( 7) was obtained by fitting the fermion gap ∆ f = 6 √ 3κ at low fields; see the inset of Fig. 2. Our mean-field results for the extended Kitaev model confirm the qualitative behavior obtained for the exactly solvable model; see Fig. 4. The charge imbalance oscillates between positive and negative values as we vary the distance from the center of the vortex, identified with the plaquette where W p < 0.Moreover, as shown in Fig. 5, the magnitude of ⟨δn l ⟩ decays exponentially with the distance from the vortex.The comparison with the result for Γ = 0 (dashed lines in Fig. 5) reveals that weak Γ and/or J interactions have an only minor effect on the ideal charge imbalance profile found in the pure Kitaev limit [30].However, changing the magnetic field direction away from the ĉ direction can induce more pronounced charge oscillations, cf.Fig. 5(b), and thus has a more substantial effect.The value of |δn l | on sites around the vortex is of the order of 10 −6 , producing local electric fields near the detection limit of state-of-the-art atomic force microscopy [30,[61][62][63].Importantly, here we use estimates for the hopping and interaction parameters for bulk α-RuCl 3 , but the charge fluctuations can be greatly enhanced if the on-site repulsion U is screened in a monolayer by the interaction with a substrate.
Since the mean-field approach allows us to treat the Zeeman term nonperturbatively, we can go beyond the results of Ref. [30] and analyze the dependence of the charge redistribution on the field direction.For a field along the ĉ direction, the charge imbalance profile is isotropic around the position of the vortex, up to small variations due to the finite distance between vortices in the finite-size system.As we tilt the magnetic field on the ac plane (perpendicular to the z bonds), a small anisotropy develops in a way that the charge imbalance is enhanced in the direction perpendicular to the field.This effect can be seen in Fig. 5(b) as the difference between ⟨δn l ⟩ for the sites that belong to the hexagon that contains the vortex (three blue dots in the center, cf.Fig. 4).
We next quantify the anisotropy in the charge distribution by computing the electric multipole moments.We note that the electric quadrupole moment has also been studied in the context of the spin nematic transition in the vortex-free ground state of a perturbed Kitaev model [64].In the limit of very large distance between vortices, the electric dipole moment vanishes because the system is invariant under spatial inversion about the vortex center.The first nontrivial multipole moment is the traceless quadrupole tensor, with components Here α, β ∈ {1, 2, 3} and R l = R l1 â+R l2 b (with R l3 = 0) is the position of site l, setting the lattice spacing to unity.Due to the finite system size, we calculate the quadrupole moment by summing over all sites within a finite radius around the vortex.This radius is taken to be slightly smaller than half the distance between vortices, but due to the exponential decay of ⟨δn l ⟩ with the distance from the vortex center, changing this radius causes only exponentially small changes in the quadrupole tensor.For a magnetic field along the ĉ direction, the rotational symmetry implies that the quadrupole tensor is diagonal and As we vary the field direction, the anisotropy is manifested in the difference between Q 11 and Q 22 and in the off-diagonal element Q 12 .Note that Q 13 vanishes identically because R l3 = 0.In a first approximation, let us discuss the dependence of the quadrupole tensor on the magnetic field direction by treating the field perturbatively in the pure Kitaev model.For magnetic field directions not perpendicular to the lattice plane, the (often discarded) contribution from second-order perturbation theory generates an anisotropic renormalization of the Kitaev couplings.This effect is captured by the Hamiltonian in Eq. ( 6) with In Fig. 6, we show the angular dependence of the quadrupole components Q 33 , Q 11 − Q 22 and Q 12 calculated from the spin correlations for the model in Eq. ( 6).
The component Q 33 does not change sign, but varies slightly around an average value with an angular dependence qualitatively similar to |h x h y h z |.In particular, Q 33 is maximum for a field along the ĉ direction, which may be interesting to maximize the intrinsic electric field produced at positions right above the vortex.On the other hand, the difference Q 11 − Q 22 vanishes for h ∥ ĉ, but is maximum when the field points along the ẑ axis; this is the direction in which the anisotropy in the effective Kitaev couplings is maximized, with Finally, Q 12 vanishes if we tilt the field along the highsymmetry ac plane, but becomes nonzero for more general field directions.
The spin correlations calculated within the mean-field approach for the extended Kitaev model lead to the same qualitative dependence on the field direction as in Fig. 6.To maximize the anisotropy in the quadrupole tensor, we focus on the direction h = hẑ, in which case all offdiagonal components vanish, and analyze how the diagonal components vary with the strength of the magnetic field.Here it is convenient to introduce the dimensionless anisotropy parameter ∆Q = (Q 11 − Q 22 )/|Q 33 |.As shown in Fig. 7, ∆Q increases with h, and the effect is more pronounced in the presence of the Γ interaction.We have also studied the case Γ < 0 and find qualitatively similar results [59].
The spatial anisotropy of the charge density profile affects the electric quadrupole interaction between vortices.Suppose the first vortex is located at the origin and the second one at r = x 1 â + x 2 b, with r = |r| much larger than the length scale in the decay of δn l .The interaction is given by the energy E of the quadrupole tensor Q (2)  of the second vortex in the electrostatic potential V (1)  generated by the first vortex, where αβ .
Since wellseparated vortices generate the same charge distribution, we now assume Q (1) = Q (2) = Q.As a result, the quadrupolar interaction can be written as When the magnetic field varies along the ac plane, the quadrupole tensor is diagonal and we obtain E = Here θ is the angle between r and â, and .Quadrupole components as a function of magnetic field direction, calculated using the exactly solvable Hamiltonian in Eq. ( 6), i.e., for J = Γ = 0.The coupling constants Kγ and κ were calculated using Eqs.( 7) and ( 9) with |h| = 0.2|K| and ∆2v = 0.263|K|.The scale is in units of t 2 2 t ′ 2 /U 3 and we set the lattice spacing to unity.Here we use L = 42.we use with the property F (−∆Q, θ) = F (∆Q, π/2 − θ).In particular, ∆Q = 0 for a magnetic field along the ĉ direction; in this case, the quadrupolar interaction becomes strictly repulsive and independent of θ.However, as illustrated in Fig. 8, the interaction can change sign for some particular directions of r if the anisotropy is strong enough.The attractive regime appears for |∆Q| > 9/7 ≈ 1.13.
According to the result in Fig. 7, this regime becomes accessible for sufficiently large h and Γ with h along the ẑ direction.We note that already in the pure Kitaev model, vortices have an effective interaction that depends on the vortex separation [65].The charge redistribution discussed here provides a mechanism to make this interaction spatially anisotropic.In the extended Kitaev model, where vortices acquire a small mobility [58], the charge density profile must be carried along with the slow vortex motion, and the anisotropic interaction may cause some nontrivial dynamics in a system of dilute vortices.Importantly, the quadrupole interaction decays algebraically with the distance between vortices; thus, at large distances it dominates over other sources of vortex-vortex interactions that are expected to decay exponentially [65].

Electrical manipulation of vortices
We now consider the effect of a local electrostatic potential on vortices.Going back to the Hubbard-Kanamori model, we couple the hole density to a potential V 0 on the six sites surrounding a hexagonal plaquette p where a vortex is located.This local potential can be generated by the electric field of a scanning tunneling microscope (STM) tip.Redoing the derivation of the effective spin Hamiltonian by second-order perturbation theory, we find that the local potential modifies the couplings on bonds between sites in p and their nearest neighbors outside p; see Eq. ( 23).In addition, the local electric potential breaks inversion symmetry and generates a Dzyaloshinskii-Moriya (DM) interaction [34].Microscopically, the DM interaction stems from crystal field splittings in the atomic Hamiltonian and asymmetries in the hopping matrix due to lattice distortions [5].We investigate this effect phenomenologically by adding to the effective spin Hamiltonian (1) the term where (αβγ) is a cyclic permutation of (xyz).The coupling D ij = D(V 0 ) is taken to be independent of the bond type γ but restricted to the bonds exterior to the plaquette with the local potential.For the DM coupling, we assume [59] such that D(V 0 ) ∝ V 0 with a dimensionless free parameter D 1 .In fact, for V 0 = 0, the DM coupling is absent since it will be generated by the tip potential.
In the solvable Kitaev model, the local electric potential lowers the energy of an isolated vortex with respect to the vortex-free configuration, but never closes the vortex gap in the absence of the DM interaction [30].In that case, this effect can be used to attract and bind vortices that have been created by some other mechanism, such as thermal fluctuations, but it does not induce vortices in the ground state of the system.Using the mean-field approach, we can now analyze how the vortex gap varies with the electric potential in the extended Kitaev model.We consider again the configuration with four equally spaced vortices, see the inset in Fig. 5(b), and apply the electric potential on the four corresponding plaquettes.The difference between the energy E 4v of this four-vortex configuration and the energy E 0v of the vortex-free state is equal to four times the vison gap.As shown in Fig. 9, the vison gap monotonically decreases with the applied electric potential, and it is further reduced for nonzero Γ and finite DM coupling D 1 .When the gap becomes too small, we encounter difficulties in the convergence of the mean-field equations.However, the extrapolation of the results indicates that the gap vanishes for sufficiently large V 0 .As a consequence, we predict that it is possible to create (or remove) vortices by modulating the local interactions, in agreement with the results of Ref. [34].We emphasize that this remarkable functionality arises due to the interplay between Γ interactions and the local DM terms induced by an STM tip.From Fig. 9, we observe that ξ 1 ∼ 0.5 is sufficient to create vortices.Using the parameters listed in Fig. 4, we find that this corresponds to realistic tip voltages of the order of V 0 < 1 V.

Extended Kitaev model
The JKΓ model in Eq. (1) follows by projecting the three-orbital Hubbard-Kanamori Hamiltonian on the honeycomb lattice, H HK = V +H so +T , to the low-energy sector spanned by a single hole per site.On-site interactions are encoded by where U is the repulsive interaction strength, J H is Hund's coupling, and the operators Ni , S i and L i are the total number, spin and orbital angular momentum of holes at site i.The operator h † iασ creates a hole at site i with spin σ ∈ {↑, ↓} and orbital α ∈ {x, y, z} for yz, xz, and xy orbitals, respectively.Defining the spinor where σ is the vector of Pauli matrices acting in spin space and l = (l x , l y , l z ) is a vector of 3 × 3 matrices that represent the effective l = 1 angular momentum of the t 2g states [30].The spin-orbit coupling term H so = λ iα h † i (σ α ⊗ l α )h i splits the degeneracy of the t 2g manifold.At each site, the low-energy subspace is spanned by the states which are associated with total angular momentum j eff = 1 2 .Finally, the hopping term in H K has the form The hopping matrix T ij in orbital space depends on the orientation of the bond between sites i and j.We label the bonds on the honeycomb lattice by γ ∈ {x, y, z} ≡ {1, 2, 3} corresponding to nearest-neighbor vectors 3 b, and δ z = − 1 √ 3 b, respectively.We parametrize the hopping matrix for a nearest-neighbor z bond as [35]  .The hopping matrix for x and y bonds follows by cyclic permutation of the orbital indices.Microscopically, the hopping parameters are associated with direct hopping between d orbitals or hoppings mediated by the ligand ions.Neglecting trigonal distortions for simplicity, hereafter we set t 4 = 0 [5,35].
The effective spin Hamiltonian for the Mott insulating phase can now be derived by applying perturbation theory in the regime U, J H ≫ λ ≫ t 1 , t 2 , t 3 .We use the canonical transformation The anti-Hermitian operator S = ∞ k=1 S k is chosen so that S k eliminates the terms that change the hole occupation numbers Ni at k-th order in the hopping parameters.We can write S k = S + k − S − k , where S + k creates excitations with Ni ̸ = 1 and S − k = (S + k ) † .For the calculation of the effective spin Hamiltonian, it suffices to consider the first-order term S 1 = S + 1 − S − 1 , with Here P (1) j is a projector onto the subspace of a single hole at site j and P (2) j,ℓ projects onto the subspace of two holes with total angular momentum ℓ ∈ {0, 1, 2}.The excited states have energies ∆E ℓ given by with J H < U/3 in the Mott insulating phase.We then take H = P low HHK P low , where is the projector onto the lowenergy subspace restricted to j eff = 1 2 states at every site.We thereby arrive at the JKΓ model [35], with an implicit sum over bond type γ, and α, β chosen so that (αβγ) is a cyclic permutation of (xyz).The couplings are In the limit t 1 , t 3 → 0 and t 2 ̸ = 0, Eq. ( 21) reduces to the exactly solvable Kitaev model [1] with a ferromagnetic Kitaev interaction (K < 0).Finally, in the presence of a potential V 0 , the respective couplings are renormalized according to , where ξ ℓ = eV 0 /∆E ℓ .

Mean-field Hamiltonian
Using the mean-field parameters in Eq. ( 5), the Majorana mean-field Hamiltonian for Eq. ( 4) is given by The first term on the right-hand side couples Majorana fermions on nearest-neighbor bonds ⟨ij⟩ γ via the 4 × 4 bond-dependent matrix The on-site term involves the matrix where V i denotes the set of nearest neighbors of site i.Finally, the constant term is We diagonalize Eq. ( 24) for N unit cells of the honeycomb lattice with periodic boundary conditions by using where c is a vector defined from 8N Majorana fermions, U is a unitary transformation, and d is a 4N -component vector of annihilation operators of complex fermions.The columns of U <(>) correspond to the eigenvectors of the mean-field Hamiltonian with negative (positive) energy.The mean-field ground state is the state annihilated by all d operators, from which we obtain the selfconsistency conditions where I = (i, µ) and J = (j, ν) combine site and fermion flavor indices.We obtain the mean-field parameters in Eq. ( 5) by setting i and j to be either nearest neighbors or the same site.Together with the mean-field Hamiltonian, Eq. ( 29) defines a set of self-consistent equations which we then solve numerically.
In our approach, we require that the constraint in Eq. ( 3) is satisfied by the mean-field solution as accurately as possible.Since ic T i G γ c i are linear combinations of operators with eigenvalues ±1, we define the quantities for the mean-field ground state average, with 0 ≤ G γ i ≤ 1.For zero magnetic field and in the absence of magnetic order, the constraints are automatically satisfied, G γ i = 0, since V 0γ i = V αβ i = 0. To describe the Kitaev spin liquid phase at finite magnetic field, we tune the Lagrange multipliers λ γ i contained in B i in order to minimize the violation of the constraint measured by G γ i .For all results shown below, we guarantee G γ i < 0.05 for all values of (γ, i).In the homogeneous case (cf.Figs. 1, 2 and 3), the largest violations occur in the vicinity of phase transitions.Away from transitions, we instead find G γ i < 10 −3 .Similarly, in the presence of vortices, the largest violations occur near a vortex but they are always bounded as specified.

Charge density coefficients
Consider the hole density operator Nl at site l in the Hubbard-Kanamori model.Using the canonical transformation in Eq. ( 18), we can write the effective charge / q t J 1 S a x / L B j B P 0 I z q Q P O S M G i v V R 7 1 S 2 a 2 4 M 5 B l 4 u W k D D l q v d J X t x + z N E J p m K B a d z w 3 M X 5 G l e F M 4 K T Y T T U m l I 3 o A D u W S h q h 9 r P Z o R N y a p U + C W N l S x o y U 3 9 P Z D T S e h w F t j O i Z q g X v a n 4 n 9 d J T X j j Z 1 w m q U H J x < l a t e x i t s h a 1 _ b a s e 6 4 = " t 4 J H 8 e C f W A n F H w r E x n z V w 9 j T D a Q = " > A A A B 6 H i c d V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q W S L 1 a C j W e q 9 8 n p j R Q K g 0 8 0 x l Q P V a / v U z 8 y + s l 2 q + 5 U x 7 G i Y a Q L R b 5 i c A 6 w t n X e M g l M C 1 S Q y i T 3 N y K 2 Z h K y r T J p m B C + P o U / 0 8 6 Z d u p 2 p X W e a l x u Y w j j 4 7 Q M T p F D r p A D X S N m q i N G A L 0 g J 7 Q s 3 V n P V o v 1 u u i N W c t Z w 7 R D 1 h v n 1 H / j U 4 = < / l a t e x i t > y < l a t e x i t s h a 1 _ b a s e 6 4 = " Y N C E J U n 0 d y p J b D i l a h x i 0 We calculate δn l using perturbation theory to leading order in the hopping matrix T ij .In systems with bondinversion symmetry like the Hubbard-Kanamori model, the first non-vanishing contribution appears at third order and is associated with virtual processes in which an electron or hole moves around a triangle [30,32,33].
To obtain this leading contribution, we generalize the hopping matrix to include hopping between next-nearestneighbor sites on the honeycomb lattice.We denote by ⟨⟨ij⟩⟩ γ a second-neighbor bond perpendicular to nearestneighbor γ bonds, see Fig. 10(a).Sizeable second-and third-neighbor hopping parameters have been calculated for Kitaev materials using ab initio methods [5,50].
For simplicity, we consider only the dominant secondneighbor hopping, which on z bonds is described by the matrix Following Ref. [30], we write the effective charge imbalance operator as δn l = (jk) δn l,(jk) , where the sum over (jk) runs over pairs of sites such that jkl forms a triangle, and each triangle is counted once.These triangles contain two nearest-neighbor bonds and one next-nearestneighbor bond, see the examples in Fig. 10(b).The calculation of δn l requires the generator of the canonical transformation up to second order in the hopping matrices, S ≈ S 1 + S 2 .After the projection onto the j eff = 1 2 subspace, we write the end result in the form Note that the effective density operator involves only twospin operators because it must be invariant under time reversal.The coefficients C µνρ jkl can be calculated as explained below.We find closed-form but lengthy expressions for general values of the hopping parameters.For t 1 = t 3 = 0 and t 2 , t ′ 2 ̸ = 0, we recover the result of Ref. [30], in which the nonzero coefficients are diagonal in spin indices, e.g., C αβ0 jkl ∼ δ αβ t 2 2 t ′ 2 /U 3 .Similarly to the derivation of the effective Hamiltonian, the addition of the subleading hopping parameters t 1 and t 3 generates off-diagonal terms in δn l,(jk) which are reminiscent of the Γ interaction.Equation (31) implies that the charge density profile of a given state is determined by its spin correlations.Charge neutrality of the Mott insulator, l ⟨δn l ⟩ = 0, implies that there is no charge polarization in a homogeneous state where ⟨δn l ⟩ is uniform.This condition is indeed satisfied when we impose that the spin correlations on different bonds respect translation and rotation symmetries, which provides a nontrivial check for the coefficients C µνρ jkl .Let us outline some steps in the calculation of the coefficients C µνρ jkl in Eq. (31).At third order in the hopping term, the canonical transformation in Eq. (30) gives δn  l,(jk) can be found in Ref. [30].The last step is to project these matrices onto the low-energy subspace spanned by the states in Eq. (17).The coefficients in Eq. ( 31) are given by Tr P low δn l,(jk) P low σ µ j σ ν k σ ρ l , where σ 0 = 1.Since the charge density operator is even under time reversal, terms that act nontrivially on an odd number of spins vanish identically, with α, β, γ ∈ {1, 2, 3}.The nonzero terms can be written as in Eq. ( 31) and depend on the specific triangle.The simplest coefficients are the ones that are already present in the solvable Kitaev model [30].For instance, for the top triangle in Fig. 10(b), we obtain where η = J H /U < 1/3.Note that this term is sensitive to the sign of the second-neighbor hopping t ′ 2 .In Ref. [30], the charge imbalance was calculated assuming a positive value of t ′ 2 , but in this work we use t ′ 2 < 0 as obtained in Ref. [50] for α-RuCl 3 .As an example for a coefficient associated with off-diagonal terms in δn l , which are generated by the hoppings t 1 and t 3 , we have

DISCUSSION
We have studied how vortices in Kitaev spin liquids generate and respond to nonuniform electric fields.While Kitaev materials are Mott insulators, charge fluctuations can be generated at low energies by inhomogeneous spin correlations that carry signatures of localized excitations.To describe this effect, we started from the three-orbital Hubbard-Kanamori model for Kitaev materials.Using a canonical transformation, we obtain effective operators in the low-energy sector in terms of spin operators that act on the pseudospin-1/2 states.The effective spin Hamiltonian is the extended Kitaev model in a magnetic field, in which the exact solvability is broken by the Heisenberg and Γ interactions as well as by a Zeeman coupling to a magnetic field.The effective density operator is generated at the level of third-order perturbation theory in the hopping terms.Generalizing the results of Ref. [30], we found that the effective density operator associated with the extended Kitaev model contains all two-spin operators allowed by symmetry, including off-diagonal terms that are absent in the pure Kitaev model.
We have developed and applied a Majorana meanfield approach which allows to consider inhomogeneous parameters.While this approach is exact for the pure Kitaev model, we have demonstrated that it captures qualitative features of the Kitaev spin liquid phase in the extended JKΓ model, where additional spin interactions are present.This model is believed to describe the candidate material α-RuCl 3 .The electric charge distribution follows by computing the spin correlations around vortices in the mean-field approach.Importantly, vortices remain localized on sufficiently long time scales even in the presence of small perturbations around the Kitaev limit, as long as the system remains deep in the Kitaev spin liquid phase.We find that the charge profile de-cays with the distance from the vortex in an oscillatory fashion.
Our results allow us to calculate the intrinsic electric quadrupole tensor of a vortex which is far away from all other vortices.The anisotropy of the quadrupole tensor can here be controlled by the magnetic field, and depending on the parameter regime, the interaction between different vortices is either repulsive or attractive.The interaction is generally enhanced by the Γ interaction.
Finally, in the presence of local STM tips near vortices, we find that one can close the vortex gap by applying a local electric potential to the tips.We thus predict that one can create vortices in a Kitaev spin liquid by means of STM tips in a controlled way.Given the recent advances in STM technology [61][62][63], our work paves the way for the electrical detection and manipulation in Kitaev materials.In particular, the successful control of Ising anyons in such materials would constitute a key step toward implementing a platform for topological quantum computation.

Figure 1 .
Figure1.Dispersion relation of Majorana fermions calculated within the mean-field approach for the homogeneous system with K = −1, J = 0, Γ = 0.2, and h = 0.4ĉ, along the indicated BZ path.For comparison, the dashed lines show the dispersion in the pure Kitaev limit (Γ = J = h = 0).

Figure 2 .
Figure2.Fermion gap as a function of magnetic field for h = hĉ along the ĉ axis, with K = −1, J = 0, and for three values of Γ: Γ = −0.1 (green circles), Γ = 0 (blue triangles) and Γ = 0.1 (red diamonds).The inset shows that for weak fields the fermion gap agrees with the perturbative result to leading order in h, ∆ f ∝ h 3 .

Figure 3 .
Figure 3. Mean-field parameter U γγ γ for fixed K = −1 and J = 0 as a function of the Γ interaction and the strength of the magnetic field along two directions: (a) h ∥ ĉ, perpendicular to the honeycomb plane; (b) an in-plane field h ∥ â. White circles represent critical points where the fermion gap closes at the Γ point in the BZ.The region labeled as KSL is identified with the Kitaev spin liquid phase.

Figure 5 .
Figure 5. Magnitude of the charge imbalance as a function of the position R1 along the zigzag path represented by the black line in Fig. 4. The dots (connected by solid lines to guide the eye) correspond to the extended Kitaev model with exchange couplings Γ/|K| = 0.2 and J/|K| = −0.02.Blue and red represent positive and negative charges, respectively.We set |h| = 0.3|K| and consider two field directions: (a) h ∥ ĉ, and (b) h ∥ ẑ.For comparison, dashed lines represent the corresponding mean-field results for Γ = 0 and otherwise identical parameters.The inset in (a) shows the corresponding case with Γ/|K| = 0.35 and J/|K| = −0.05for h ∥ ĉ (filled circles), comparing with the results for Γ/|K| = 0.2 and J/|K| = −0.02 in the main plot (empty circles).The values of δn l are in units of |t 2 2 t ′ 2 /U 3 |.The inset in (b) shows the geometry with four equally spaced vortices on the torus with a smaller system size.

Figure 6
Figure 6.Quadrupole components as a function of magnetic field direction, calculated using the exactly solvable Hamiltonian in Eq. (6), i.e., for J = Γ = 0.The coupling constants Kγ and κ were calculated using Eqs.(7) and (9) with |h| = 0.2|K| and ∆2v = 0.263|K|.The scale is in units of t 2 2 t ′ 2 /U 3 and we set the lattice spacing to unity.Here we use L = 42.

Figure 9 .
Figure 9. Energy of the four-vortex state vs applied electrostatic potential V0, with the dimensionless quantity ξ1 = eV0/(U − 3JH ), for different values of Γ and of the DM coupling D1; see Eq. (14).Symbols represent mean-field results for the extended Kitaev model with K = −1, J = 0 and the magnetic field h = 0.2 The linear system is L = 28, and solid lines are a guide to the eye only.They were obtained by a fit to the function a + b tanh[c(ξ1 − d)].
k X Z u y 5 f 1 S 5 L l b s s j j y c w C m c g w c 3 U I E H q E I d G C A 8 w y u 8 O U P n x X l 3 P h a t O S e b O Y Y / c D 5 / A N R B j P g = < / l a t e x i t > j < l a t e x i t s h a 1 _ b a s e 6 4 = " d / R h c R S M b E v 4 i v p b A n d L 8 T o 6 s f g 1 g g P A M r / D m P D o v z r v z M W 9 d c f K Z I / g D 5 / M H 1 c W M + Q = = < / l a t e x i t > k < l a t e x i t s h a 1 _ b a s e 6 4 = " F k O 2 T A I B l l + W e e d 9 S p S u K P Y E V / A = " > A A A B 6 H i c d V D J S g N B E O 2 J W 4 x b 1 K O X x i B 4 G m Y m q 7 e g F 4 8 J m A W S I f R 0 a p I 2 P Q v d P W I I + Q I v H h T x 6 i d 5 8 2 / s S S K o 6 I O C x 3 t V V N X z Y s 6 k s q w P I 7 O 2 v r G 5 l d 3 O 7 e z u 7 R / k D 4 / a M k o E h R a N e C S 6 H p H A W Q g t x R S H b i y A B B 6 H j j e 5 S v 3 O H

Figure 10 .
Figure10.Honeycomb lattice with nearest-and next-nearestneighbor hopping.(a) Red, green and blue lines correspond to γ = x, y, z bonds, respectively.The hopping matrix on nearest-neighbor bonds (solid lines) is written in terms of hopping parameters t1, t2, t3, and t4.On second-neighbor bonds (dashed lines), we consider a single hopping parameter t ′ 2 .(b) Examples of triangles that contribute to the effective charge density operator at site l.

2 [
Nl , S + 1 ] + h.c., where we organize the contributions in terms of triangles with site l at one vertex.In this notation, the contribution from each triangle with two other sites (jk) ≡ (kj) contains two terms, δn(3)l,(jk) = δn . Explicit expressions for the matrix elements of δn(3)
The corresponding matrices for x and y bonds follow by cyclic permutation of the indices.Assuming |t ′ 2 | ≪ |t 1 |, |t 2 |, |t 3 |, we calculate the charge density response to first order in t ′