Anyons in integer quantum Hall magnets

Strongly correlated fractional quantum Hall liquids support fractional excitations, which can be understood in terms of adiabatic flux insertion arguments. A second route to fractionalization is through the coupling of weakly interacting electrons to topologically nontrivial backgrounds such as in polyacetylene. Here we demonstrate that electronic fractionalization combining features of both these mechanisms occurs in noncoplanar itinerant magnetic systems, where integer quantum Hall physics arises from the coupling of electrons to the magnetic background. The topologically stable magnetic vortices in such systems carry fractional (in general irrational) electronic quantum numbers and exhibit Abelian anyonic statistics. We analyze the properties of these topological defects by mapping the distortions of the magnetic texture onto effective non-Abelian vector potentials. We support our analytical results with extensive numerical calculations.

Strongly correlated fractional quantum Hall liquids support fractional excitations, which can be understood in terms of adiabatic flux insertion arguments. A second route to fractionalization is through the coupling of weakly interacting electrons to topologically nontrivial backgrounds such as in polyacetylene. Here we demonstrate that electronic fractionalization combining features of both these mechanisms occurs in noncoplanar itinerant magnetic systems, where integer quantum Hall physics arises from the coupling of electrons to the magnetic background. The topologically stable magnetic vortices in such systems carry fractional (in general irrational) electronic quantum numbers and exhibit Abelian anyonic statistics. We analyze the properties of these topological defects by mapping the distortions of the magnetic texture onto effective non-Abelian vector potentials. We support our analytical results with extensive numerical calculations.

I. INTRODUCTION
The discoveries of the integer and fractional Quantum Hall (QH) effects have revolutionized condensed matter physics: the important concept of a topological invariant was introduced to explain the quantized Hall conductivity of the former [1], while the novel notion of topological order, i.e., a type of nonlocal order with no Landau symmetry-breaking and no local order parameter, was introduced to describe the latter [2]. Topological order goes hand in hand with exotic phenomena such as fractional charge and statistics [2]. While the strongly correlated, topologically ordered fractional QH systems indeed have fractional quasiparticles [3,4], their more traditional weakly correlated counterparts have quasiparticles with integer charge.
Here we show that anomalous integer QH systems, which can emerge even in the absence of external magnetic fields in frustrated magnets, provide a new playground for electronic fractionalization. While in conventional integer QH systems (two-dimensional electron gas in a magnetic field), fractional charge can only be induced artificially, e.g., near superconducting vortices [5], we demonstrate that intrinsic topological defects in the noncoplanar magnetic systems, which exhibit anomalous integer QH effect, naturally harbor excitations with fractional electronic quantum numbers. We also show that these defects exhibit an Abelian anynonic exchange statistics.
The existence of fractional excitations in a QH liquid can be deduced from Laughlin's argument: upon adiabatic local "insertion" of a flux quantum, a fractional charge q = σ xy e, where σ xy is the fractional Hall conductivity (in units of e 2 /h), flows in from infinity [3,6]. Since integer flux quantum can be "gauged away" if its core is smaller than the physical unit cell size, the charge q is in fact the charge of an elementary quasiparticle. Naturally, insertion of a fractional flux can also lead to the same result with an integer Hall response. Fractional fluxes, however, cannot be gauged away, and thus cannot be naturally associated with quasiparticles. It may thus appear that the existence of intrinsic fractional excitations in a QH system requires a fractional σ xy .
An alternative route to electronic fractionalization is via interaction between electrons and topological defects in some order parameter, such as a pattern of lattice distortions in polyacetylene or graphene-like structures [7][8][9][10][11][12], a superconducting vortex [13][14][15][16], or a meron in a QH bilayer systems [17,18]. Here we analyze the scenario where QH effect itself emerges as a result of the coupling of itinerant electrons to chiral magnetic states. Such states can form through spontaneous time-reversal symmetry breaking [19][20][21][22][23][24][25], as well as with an assistance of an external magnetic field that explicitly breaks the time-reversal symmetry. Remarkably, the topologically stable defects in such ordered magnetic media act as effective fractional fluxes, giving rise to natural excitations with fractional charge. (In a broader context, studying topological defects in topological phases has attracted considerable recent arXiv:1306.6080v1 [cond-mat.str-el] 25 Jun 2013 interest [26][27][28][29][30].) The key ingredients of our model system are i) localized magnetic moments capable of forming a noncoplanar state that can be smoothly distorted at low energy cost, and ii) itinerant electrons which interact with the local moments, and possibly induce this state. Unlike collinear (such as the Néel) and coplanar (such as the "120-degree") order, noncoplanar magnetic states are rarely stabilized for classical magnetic moments with short-range interactions (see Ref. [31] for an exception). However, when magnetic moments are coupled to itinerant electrons, noncoplanar states are quite common, i.e., the magnetic phase diagrams of such itinerant systems seem to generically contain energetically stable phases with noncoplanar magnetic ordering [21,23,24,[32][33][34][35][36][37]. The interplay of noncoplanar magnetic moments and itinerant electrons can then lead to spontaneous quantized integer quantum Hall effect as a result of the nontrivial Berry phases imparted to the electrons by the noncoplanar magnetic texture [21,23,24]. Such spin-chirality-driven Hall effect may be realized in a wide range of materials such as mangenites, CrO 2 , the element Gd, the cobaltates, and pyrochlore ferromagnets (see Ref. [38] and the references therein). To illustrate the generic nature of our results, here we consider two different models -the triangular and kagomé Kondo lattices -both exhibiting the same kind of electronic fractionalization.
The outline of this paper is as follows. In Sec. II, we summarize our main results. In Sec. III, we briefly introduce the two models on the triangular [24] and the kagomé [21] lattice. In Sec. IV, we discuss the topologically stable defects of noncoplanar magnetic states and give a general argument for charge fractionalization. In Sec. V we present a microscopic theory of the fractional charge in limit of large Kondo coupling in the triangular lattice model, where the system maps to a simple model of spinless fermions. In Sec. VI we present numerical results for the properties of individual vortices, confirming the analytical calculations. In Sec. VII, we discuss the mutual exchange statistics of the vortices and present extensive numerical calculations, verifying our theoretical prediction of a relationship between the exchange statistics and the charge and magnetization of the vortices.

II. SUMMARY OF RESULTS
Our results concern systems where local moments S i form noncoplanar magnetic states, such as in Figs. 1(a) and 1(b). When itinerant electrons interact with these magnetic states, for certain densities and ranges of interaction strength, they exhibit anomalous QH effect. The generic Hamiltonian that leads to this behavior has the form: where electrons hop on a two-dimensional (2D) lattice and interact with the local magnetic moments via onsite exchange interaction. Summation over repeated site (roman) and spin (greek) indices is implied. Here t i j is the intersite hopping (in this work we only consider nearest-neighbor hopping t i j ), J is the exchange interaction constant, S i is the local magnetic moment (assumed classical), σ = (σ x , σ y , σ z ) is the vector of Pauli matrices, and c iα is the operator of electron annihilation on site i with spin α. H S is the classical Hamiltonian that only includes spin variables. The existence of topologically stable vortex defects in a noncoplanar magnetically ordered medium follows from the nontrivial fundamental homotopy group of the space of energetically degenerate configurations (order-parameter space). i) on the triangular lattice with H S = 0, the time-reversal symmetry can be broken spontaneously, leading to a spontaneous QH ground state (Fig. 1a.) for several electron densities [24,[32][33][34]. As the rotation of this texture around any axis keeps the energy unchanged, the the order-parameter space is the S O(3) group. In this case, the nontrivial fundamental homotopy group π 1 (O(3)) = π 1 (S O(3)) = Z 2 [39] guarantees the existence of Z 2 vortices. ii) on the kagomé lattice with H S representing Heisenberg interactions between the local moments and their coupling to an external Zeeman field, the timereversal symmetry is explicitly broken by the Zeeman field. In this model, the degenerate ground-state manifold has S O(2) symmetry (as the texture can be rotated around the axis determined by the Zeeman field), and π 1 (O(2)) = π 1 (S O(2)) = Z results in vortices characterized by an integer winding number.
An example of the vortex spin texture is shown in Fig. 1c  and 1d. For instance, it can be obtained by rotating the order parameter (every magnetic moment) by a position-dependent angle around theẑ axis (assuming the vortex core is at the origin); ν is the winding of the vortex. In the triangular (kagomé) lattice case, the two topological classes set by the parity (value) of ν. Unless stated otherwise, here we consider the ν = 1 case. Our main result is that the topologically nontrivial vortices have a fractional electric charge. If the vortex has a fixed axis of rotationn (in the kagomé lattice case above the axis is pinned to the magnetization direction), the value of the fractional charge is given by where m is an integer, and σ 00 xy is the usual QH conductance, which characterizes the charge current flowing in the x direction in response to an electric field in y direction. Similarly, σ 0n xy characterizes the charge current flowing in the x direction in response to a "spin-n electric field" (to be defined later on) in the y direction. For the textures on the triangular lattice without a uniform magnetization, the latter offdiagonal responses vanish, and for an integer QH system with σ 00 xy = e 2 /h, the charge is given by independently of the axis of rotation. We further show that vortices with a fixed axis of rotation, carry a net magnetization m n in the direction of this axis, which stems from the spins of itinerant electrons. We also demonstrate that vortices have anyonic exchange statistics, with a statistical angle that is related to both the charge and the magnetization of the defect, which for q given by Eq. (4) is where p is an integer and m n is the aforementioned magnetization of the vortex. In addition to analytical arguments, we perform extensive numerical computations to test our results.

III. MODELS AND INTEGER QUANTUM HALL EFFECT
The energy of classical configuration of magnetic moments in Eq. (1) depends on the quantum electrons. The configuration of the magnetic moments can be thought of as a set of external parameters for a quadratic electronic Hamiltonian, which uniquely determines the electronic energy at zero temperature for a given filling fraction.
On the triangular lattice with nearest neighbor hopping, it has been shown that the magnetic moments spontaneously form an "all-out" noncoplanar texture around 1/4 and 3/4 filling fraction. At 3/4 filling, this instability can be understood in terms of the perfect nesting of the Fermi surface with the three ordering wave vectors corresponding to the all-out state [24]. For general fillings the lowest energy states can be obtained by Monte-Carlo simulations [32][33][34]. At the filling fractions that correspond to electronically gapped all-out state, the system exhibits an integer quantum Hall response.
Similarly, integer QH appears on the kagomé lattice with the moments forming an umbrella-like state where the spins are canted away from a 120-degree ordered state [21], at filling fractions 1/6, 2/6, 4/6 and 5/6. This state can be either induced by an external magnetic field in magnetic systems with the nearest-neighbor antiferromagnetic exchange interaction, or can emerge spontaneously, due to spin-orbit coupling. Both ways correspond to specific forms of the classical term H S in Hamiltonian (1). In the case of nonzero external magnetic field, which we consider here, H s has a S O(2) degenerate ground-state manifold characterized by i∈ S i ∝ H, where H is the magnetic field, and i∈ indicates a vector sum of the magnetic moments in a triangle. Classically, order by disorder selects coplanar states at infinitesimal temperatures [40]. By generating a large number of such coplanar states satisfying the above constraint, and computing the ground-state energy of the fermionic Hamiltonian, we verified that these generic configurations have a higher energy than the "umbrella"-like noncoplanar state of Fig. 1 at zero temperature. This indicates that in the presence of itinerant electrons, this noncoplanar state is selected out of the classical degenerate manifold. In the following we therefore assume that the spins form an umbrella state with the canting angle arctan 1 2 √ 2 with respect to the plane of the lattice (chosen so that, as shown in Figs. 1a and 1b, the three magnetic moments point to three corners of a tetrahedron similarly to the all-out structure of the triangular lattice case).
The integer quantum Hall response in both cases described above stems from interaction between electrons and the noncoplanar magnetic state. Electrons hopping in a noncoplanar magnetic state are subject to an effective Berry phase, which results in a gapped integer quantum Hall liquid. This can be understood most simply in the limit of large exchange coupling J [41]. In this limit, the problem can be projected to a spinless hopping model, with each spinless electron representing an electron that is spin-polarized in the direction of the local magnetic moment. The (gauge-dependent) hopping amplitude can be constructed as shown in Fig. 2. Even though each individual hopping phase is gauge-dependent, the flux Φ 123 = arg χ 1 |χ 2 χ 2 |χ 3 χ 3 |χ 1 through every triangular plaquette is gauge-invariant, and is given by half the solid angle subtended by the three magnetic moments S i [42]. This flux is generically nontrivial for noncoplanar textures: in the particular case of Ref. 24 (see Fig. 1), e.g., we have a flux π/2 through each triangular plaquette. This Berry phase has a similar effect to an external magnetic field and gives rise to bands with nontrivial Chern numbers. As seen in Fig. 3, the spectrum of the triangular lattice model consists of four doubly degenerate bands with Chern number +1, −1, −1, and +1. Similarly, as seen in Fig. 3, the spectrum of the kagomé lattice model [21] consists of six bands with Chern number −1, 0, +1, +1, 0 and −1.

IV. TOPOLOGICAL DEFECTS AND FRACTIONALIZATION
A. Topological defects from homotopy theory As mentioned before, the spin-rotational symmetry can be broken either explicitly or spontaneously, which affects the structure of the degenerate manifold ("the order parameter space"). In our kagomé lattice model, the out-of-plane direction of the magnetic field is fixed, but all the magnetic moments can be simultaneously rotated around this axis while preserving energy. Therefore the order parameter space is given by S O(2), which can be geometrically represented by a unit circle. A loop in real space then maps to a loop in the order-parameter space, giving rise to regular S O(2) vortices, which as illustrated in Fig. 4(a), are characterized by an integer winding number π 1 (S O(2)) = Z.
On the other hand, in the triangular lattice model there is no such preferred axis in the absence of an external Zeeman field or spin-orbit coupling. Then the order-parameter space corresponds to a full rotation matrix in 3D, which can be parametrized by an angle and an axis. This space can be geometrically represented by a solid sphere of radius π with antipodal points on the surface identified: the distance from each point to the center of the sphere represents the angle of rotation, while the vector connecting the point to the center gives the axis of rotation. The identification of antipodal points follows from the fact that clockwise and counterclockwise rotations around the same axis by angle π are equivalent. A 1D loop in real space maps onto a 1D loop in the order-parameter space, which as seen in Fig. 4(b) can fall into two distinct topological categories: contractible (topologically trivial), and noncontractible (topologically nontrivial). Mathematically, this classification is encoded in the fundamental homotopy group of the 3D rotations π 1 (S O(3)) = Z 2 [43].
The noncontractible loop corresponds to a nontrivial vor- tex. For example, the vortex configuration of Fig. 1(d), which has a well-defined axis of rotation in theẑ direction, corresponds to the following noncontractible loop: a straight line passing through the North pole, the center of the sphere and the South pole (which is identified with the North pole) as shown in Fig. 4(c).

B. Fractionalization from Laughlin's argument
Vortices are inhomogeneities in the magnetic texture, which correspond to a position-and possibly time-dependent distortion of some reference state, S i = R(r i , t)S 0 i . These inhomogeneous states can be mapped onto a state with homogeneous order parameter, but in the presence of an effective (in general) non-Abelian vector potential as follows. The rotation of the order parameter in the classical spin space can be transformed into a unitary rotation U(r i , t) ≡ U i of the electron spinors, This mapping allows one to conveniently calculate the charge and spin currents in response to the order parameter distortions. In particular, the vortex configuration corresponds to a spatially localized non-Abelian flux, which can be used to determine the charge of the vortex.
Assuming that the variation of the texture is slow on the lattice constant scale, we can make an expansion, a σ a , with the indices ν = {t, x, y} representing the space-time components and a = {1, 2, 3} the S U(2) generators. The Hamiltonian can then be written as where H 0 is the Hamiltonian corresponding to the static undistorted spin state, i.e., Eq. (1) with the substitutions c → ψ and S i → S 0 i , and the currents are defined as Denoting the 2 × 2 unit matrix by σ 0 , for a = 0 the definitions above also give the charge density and current operators [44].
What is the vector potential that corresponds to a vortex? The vortex texture is necessarily singular near the core; therefore, the transformation to unwind the vortex is singular as well. The simplest transformation that takes the vortex texture (with the z axis as the axis of rotation) into a uniform one is e iσ 3 φ(r)/2 , where φ(r) is the angle of rotation [see Eq. (2)] around the z axis. However, since upon going around the vortex, φ(r) → φ(r) + 2π, the unitary changes sign, this would correspond to antiperiodic boundary conditions for fermions along a line connecting the vortex to infinity. To avoid this complication, the above S U(2) transformation can be augmented by a U(1) one [45]. The combined transformation U(r) = e i(σ 3 +1)φ(r)/2 is only acting on up-fermions. Its associated vector potential is which has the field strength zero everywhere except for the vortex core. Due to the singular nature of the vector potential, we cannot directly apply the linear-response formalism in the vicinity of the vortex core. To calculate the vortex quantum numbers, we can instead invoke an analog of the Laughlin's argument [3,6]. The flux of the non-Abelian gauge-like field through the vortex core is Now, suppose that the flux is turned on adiabatically from zero to Φ. That will generate a non-Abelian emf acting on electrons, which at large distances from the core will be nearly uniform (tangential to any circle centered at the vortex). The vortex quantum numbers are then obtained by integrating the associated currents generated in response to this emf. If the texture is slowly varying in time (compared with the inverse energy gap in the spectrum), as well as in space, the vector potential A ν a is small and the expectation values of the current operators defined in Eqs. (9) and (10) can be calculated with the linear-response theory. The charge and spin current are related to the vector potential through J η a = σ ab xy ηµν ∂ µ A ν b . The Laughlin argument yields the charge by integrating the current in a dynamical process where the vortex is created adiabatically (the lattice provides an underlying regularization). This argument relies on two conditions: First, we need to have a continuous sequence of gapped Hamiltonians connecting the one with flux zero to the one with Φ. Second, we need a continuity equation relating the currents we can calculate in linear response to quantum number densities. We have explicitly identified a sequence of gapped Hamiltonians in the limit of large J and we thus expect that an adiabatic process exists for an arbitrary J as well. The quantum numbers of interest for us are charge and magnetization. Since total electron number commutes with the Hamiltonian, the charge current strictly satisfies a continuity equation and we can use the Laughlin argument to compute the charge quantum number. We do not have such continuity equation for the spin current, and thus our vortex defects do not have a well-defined spin quantum number that would be independent of the exact vortex configuration.
We are now in a position to state our main result [Eq. (3)] for the fractional charge. If we have a fixed axis of rotation, sayẑ as in Eq. (11), the charge current flowing toward the vortex core gets contributions from σ 00 xy = −σ 00 yx and similarly from σ 03 xy and σ 03 yx . Note that in general the underlying lattice could break rotation symmetry and the last two response functions could be different. Using Laughlin's argument, we immediately obtain the second term in Eq. (3): An "electric field" E in the tangential direction gives a current J = (σ xy E cos θ, −σ yx E sin θ) at polar angle θ. The current flowing toward the vortex core is then given by J r = −σ xy E cos 2 θ + σ yx E sin 2 θ. Since the average of both sin 2 θ and cos 2 θ is 1/2, Eq. (3) (modulo the integer m) follows upon integration over θ.
The origin of an undetermined additive integer can be understood as follows. There are many possible choices of the single-valued gauge transformations that unwind the vortex; e.g., another choice could lead to the vector potential A ν = (−1 + σ 3 ) ∂ ν φ/2. While such different choice does not change the magnetization, the charge accumulated in the vortex core changes sign. This ambiguity is naturally understood in terms of the electron occupancy of a particular localized electronic state, ε 0 , inside the spectral gap. When this state is empty, the charge of the system is q, and when it is occupied, the charge is q + 1, all relative to the uniform state. In general, there can be more than one localized state inside the vortex core. Occupying any of these states increases the vortex charge by one electron charge [this corresponds to more general choices, A ν = (1 + 2n + σ 3 ) ∂ ν φ/2, with n any integer]. In case when the order parameter space is S O(3), as is the case for triangular lattice model in zero magnetic field, by direct calculation we can verify that σ 0a yx vanishes for a 0 and the charge of the vortex remains half-odd-integer for an odd vorticity. On the other hand, for an even vorticity, the charge induced according to the Laughlin argument will be integer. For Z 2 vortices, this is consistent with the homotopy classification that says that double vortex can be smoothly connected to a uniform state, and thus the charge of quasiparticle associated with the double vortex can only be integer.
Through a explicit calculation described in Appendix. B, we can compute the necessary linear-response functions. As stated before, for the triangular lattice model we find: triangular : σ 00 xy = −σ 00 yx = −e 2 /h, σ 0a xy = σ 0a yx = 0, a 0 (12) The signs of both conductivities are flipped by switching between 3 / 4 and 1 / 4 fillings, or by changing the sign of the chiral ordering. The kagomé lattice model, on the other hand, has net magnetization, and, hence, in addition to σ 00 xy = −σ 00 yx = e 2 /h, the following off-diagonal responses are nonzero: This results in the same q = 1/2 charge for an axis of rotation in the xy plane (σ 01 xy = σ 10 xy ≈ 0), but if the axis of rotation has a component in theẑ direction (the direction of the overall magnetization), we get nonuniversal J-dependent contributions to the fractional charge.
Additionally, the adiabatic creation of the vortex generates a spin current through nonvanishing responses such as σ aa xy for a 0. For example, in the triangular lattice model we find σ aa xy = 1 3 e 2 h . As we stated before, there is no continuity equation for spin. However, since the divergence of the induced current is zero far from the vortex core, it is expected that the spin current will be nearly conserved everywhere, except near the vortex core where spin density accumulates. This suggests that the magnetization attached to a vortex might be still close to the expected value obtained from integrating the spin current. With the assumptions above, since the final flux both in spin-σ 3 and charge channels is half of the flux quantum, the accumulated magnetization for the triangular lattice may be close to [46] m z ≈ 1/12. (the extra 1/2 for m z ≈ 1 2 × 1 2 × 1 3 is due to the fact that electron spin is σ/2.) We will numerically examine this approximate result for m z in the subsequent sections, and find that it underestimates the average magnetization by up to 40% due to spin nonconservation.

V. MICROSCOPIC DERIVATION FOR J → ∞
In this section we present a direct microscopic derivation of the fractional charge in the limit of J → ∞. We only consider the triangular lattice model for brevity. This calculation is illuminating as it provides a step-by-step derivation of the charge accumulation directly on the lattice. As stated before, due to the alignment of the electron spin with local moments, in this limit we have a model of spinless fermions coupled to U(1) fluxes (for arbitrary values of J we had spinful fermions coupled to an S U(2) fluxes). Therefore, this calculation also sheds light on the mechanism for the emergence of fractional flux.
The t i j → t i j = t χ i |χ j mapping, where |χ i is a spinor in the S i direction, modifies both the magnitude and phase of the hopping amplitude. As we will see, however, the important physics stems from the change in phase. The following effective Hamiltonian then captures the physics of the problem in the large-J limit: We pick a coordinate system to represent the magnetic moments of Fig. 1(a) (this choice is arbitrary as there is no spinorbit coupling). Consider a chiral texture of magnetic moments on the triangular lattice as in Fig. 1(a) [we have written out explicit spin components in Fig. 5(a)   In the basis Ψ † k = c † 1k c † 2k , with the annihilation operators c 1 and c 2 shown in Fig. 5(b), the Hamiltonian can be written as a 2 × 2 matrix H 0 (k) in momentum space, i.e., obtained by rotating each magnetic moment S i around theẑ axis by an angle θ equal to the polar coordinate of site i (in a planar polar coordinate system with the vortex core at the origin). The rotation angles are not small but their difference for two nearby magnetic moments is small far away from the vortex: it scales as r −1 . We have chosen a convenient gauge here such that the spin rotations due to the vortex give a small perturbation to the Hamiltonian far away from the vortex core.
Our approach is then to calculate the expectation value of the current flowing toward the vortex core, in a region of size depicted in Fig. 6, which is much larger than the lattice spacing a and much smaller than its distance r from the vortex core, while turning on the vortex adiabatically. Note that the lattice provides short-distance regularizations so the singularity of the vortex at the core is not an issue (unless the core is sitting right on a lattice site). The condition a r allows us to treat the problem in the region of size as approximately translationally invariant (analogous to gradient expansion methods [10,47]). Through a Laughlin-type argument, we can then find the charge bound to the vortex by integrating this current.
The first step to carrying out the above procedure is to find the first-order correction to φ ab in the region shown in Fig. 6. The correction for each bond depends on the difference between the rotation angles of the magnetic moments at the two ends of the bond, which we represent by δ ab . To leading order, in a region labeled by r and θ as in Fig. 6, δ ab depends only on the direction of the bond. There are three types of bonds corresponding to the three lattice vectors a i , i = 1, 2, 3 so we get three types of δ ab ≡ δ i , i = 1, 2, 3: We can then compute the correction to φ ab to first order in δ i using simple Taylor expansions. The results of these calculations are represented graphically in Fig. 5(c). We observe that even in the presence of the variations δ i , the Hamiltonian can be written as a 2 × 2 matrix in the same basis as before: where V(k) is the perturbation to the Hamiltonian due to the presence of the vortex (which depends on θ and r and is valid to leading order in the region shown in Fig. 6).
Moreover, we can represent the current operators J x (k) = ∂H 0 (k) ∂k x and J y (k) = ∂H 0 (k) ∂k y , as 2 × 2 matrices and write the tangential and radial current operators as J θ = −J x sin θ + J y cos θ, J r = J x cos θ + J y sin θ.
We are now ready to state the key result of this section. By explicitly writing out V(k) (see Appendix A for details), we find a relationship between the current operator and the perturbation to the Hamiltonian: This relationship was derived in a fixed (global) gauge and at an operator level. Let us now consider a flux π ≡ −π inserted locally in the system. In the continuum limit, this local flux corresponds to a tangential vector potential A θ = π/2πr = 1/2r along a circle of radius r and perimeter 2πr in some gauge. As A θ also couples to J θ , we find that the vortex effectively acts like a (fractional) flux π. Now, we know that if a flux π in adiabatically inserted into an integer quantum Hall system with σ xy = e 2 /h, a fractional charge 1/2 will be transported to the flux tube according to the Laughlin's argument.
Since the charge is a half the Hall conductance -which is a topological invariant -any small variation in the magnitude of the hopping amplitude, which we neglected earlier, can not change it. One subtle issue with the above argument is that although we have a (gauge-dependent) operator relationship suggesting that the vortex effectively acts as a fractional flux, the gaugeinvariant fluxes induced by the vortex are completely different from a localized flux π. In fact the vortex corresponds to an intricate pattern of fluxes that only decays with the distance from the vortex core as 1/r. The fractional charge does not correspond to a well-defined magnetic flux bound to the vortex. This means that the total flux through any closed loop around the vortex depends on the geometry of the loop and does not converge by increasing the loop size. Despite this intricate flux pattern, we can perform a direct linear-response calculation for a dynamical process where the perturbation V(k) is turned on adiabatically, and as expected, we do obtain q = 1/2. This direct calculation is presented in Appendix A.

VI. NUMERICAL RESULTS FOR CHARGE AND MAGNETIZATION
We now check the above results numerically. We plot, in Fig. 7, the charge and magnetization distribution in the vicinity of the vortex core for the triangular lattice model (the kagomé lattice gives similar distribution). As expected, the charge localized in the core is half-odd-integer for odd winding and integer for even winding. The agreement between the vortex magnetization obtained numerically with what the Laughlin's argument would give is not perfect because the spin current is not conserved. Nevertheless, particularly for large J, the discrepancy is not too large and we verify that the vortex spin polarization approximately scales with the vortex winding, as shown in Fig. 8(b).
We have also considered deviations from the fully symmetric assumptions. In particular, we added a Zeeman field h alongẑ axis acting on electrons, i.e., H → H + h i c † iα σ z αβ c iβ . We verified that as long as the spectral gap does not close, the charge Hall conductivity does not change. However, in contrast to the fully symmetric case, new nonzero response functions emerge, namely σ 03 xy , σ 03 yx 0, which correspond to the charge response to A 3 , the σ 3 component of the vector potential. In Fig. 8(a), we plot the dependence of the vortex charge on h. The solid line is q = −σ xy /2 − (σ 03 xy − σ 03 yx )/4 − 1, which directly follows from the application of the Laughlin argument. Note that the charge is generally irrational and determined modulo an integer. The agreement is very good, all the way to the value of h where the gap in the electronic spectrum closes.  For the kagomé lattice case, we numerically computed the charge bound to a vortex withn =ẑ rotation axis, which is the only allowed axis for a bare texture with magnetization in the z direction. Fig. 9(a) shows the charge of a vortex in the kagomé lattice for several coupling strengths. The horizontal axis does not have a linear scale so both the variations at small J and the saturation for large J can be displayed. The results computed from the Laughlin argument show very good agreement with the numerical ones. Interestingly, putting an ad hoc vortex with thex axis of rotation also gives a charge consistent with the Laughlin argument (such ad hoc vortex will not be stable and the charge will change once the vortex relaxes.) Fig. 9

VII. EXCHANGE STATISTICS
We now turn to the exchange statistics of the vortices, considering the case of triangular lattice in zero Zeeman field for concreteness. There are two distinct possibilities for a combination of two vortices: (1) the total charge of two vortices is even, or (2) the total charge is odd. The former case will be realized if the vortices are pulled apart from the uniform "vacuum": since the initial state has total charge zero (relative to the uniform background), the state with two vortices will keep the same charge. Since the charge of an individual vortex is half-odd-integer, for large inter-vortex separation there are two energetically equivalent ground states that correspond to vortex charge configurations ( 1 / 2 , − 1 / 2 ) and (− 1 / 2 , 1 / 2 ). This degeneracy can lead to non-Abelian exchange statistics, but unlike the Majorana states in the p-wave superconductorlike systems [15,16], there is no topological protection. In other words, any local disorder can shift the bound state energy in a given vortex and lift the degeneracy.
The case (2) can be obtained, e.g., upon electron or hole doping of vortex bound states in the system with equally charged (half-odd-integer) vortices. In this case, there is no ground-state degeneracy and exchanging two vortices can only give rise to an Abelian Berry phase. We argue that in this case (absence of degeneracies), the vortices have anyonic statistics with a phase determined by the fractional charge and the magnetization of the vortex. Since magnetization can depend on microscopic details such as the location of the vortex core with respect to the lattice, the statistical angle has pathdependent contributions. Despite these subtleties, we verify in this section, through extensive numerical calculations, that the statistical angle is indeed linearly related to the magnetization.
Consider vortices obtained by spin rotations along theẑ axis. As argued in the previous sections such vortices produce a non-Abelian flux Φ = Adr = (1 + 2m + σ 3 )π, which is a 2 × 2 diagonal matrix coupled to spin-dependent density n = n ↑ n ↓ . Therefore we expect the total Berry phase accumulated upon adiabatically taking one vortex all the way around the other to be equal to Noting that q is half-odd integer, we obtain tr (Φn) = pπ + π/2 + 2πm z , where p is an integer, which depends on microscopic details such as the occupation of midgap states. Let us comment that the above relationship is an exotic feature of noncoplanar textures coupled to spinful electrons. In the more traditional case, often discussed in the literature, fractional particles of charge 1/2 usually have statistical angle π/4. Indeed from our discussion in Sec. V, we may naively expect such statistics. In the spinless model, we found that one vortex effectively acts as a U(1) flux π. Now taking another vortex of charge 1/2 around this vortex should result in a Berry phase of π/2 and Θ = π/4. What is wrong with this argument? In the spinless case, we write a model and construct wave functions in a basis labeled by lattice sites. Each of the lattice sites, however, represent the projection of the electron spin along the direction of the magnetic moment. As long as the vortices are static, we do not need to worry about this additional information, which is lost in the spinless model. However, if we have moving vortices, working with a spinless model amounts to working in a time-dependent basis.
To verify Eq. (16), we numerically compute the statistical angle through exact diagonalization (for technical details see Appendix C). To find the statistical angle, we compute the Berry phase φ bare acquired by taking a vortex along a closed path, which does not enclose another vortex, and the Berry phase φ vortex obtained by taking it around another vortex (on the same path). The statistical angle is then given by [12] θ = (φ vortex − φ bare )/2.
For efficiency, we have done our numerics in the limit of large J. This allows us to diagonalize smaller matrices (by a factor of 2), but to compute the overlaps we need to put back in the information regarding the local magnetic moments and expand our spinless wave functions in a spinful basis: if we have amplitude ϕ i on site i in the spinless wave function, and local moment S i on that site, the amplitude in the spinful basis is simply ϕ i |χ i , where |χ i is a spinor in the direction of S i .
We have performed our calculations for three system sizes L = 40, 50, 60 lattice spacing (see Fig. 12 ) with open boundary conditions (for brevity we only present data for L = 60 but similar conclusions can be drawn from smaller sizes). We have used path radii r/L = 0.22, 0.24, 0.26, 0.28, and for each system size and path radius, we have used several numbers of particles all inside the quarter-filling spectral gap with unoccupied and occupied vortex-bound midgap states. We found that, as expected, changing the particle number by filling edge modes does not affect the Berry phase. Thus we present results for only two particle numbers: one with all vortex bound states empty and one with a filled bound state. We have performed our calculations with 512 and 1024 discretization point, and obtained good convergence. The results are shown in Table. I.
The results above do show path-dependent fluctuations of around 10−15%. However, there is good stability for different system sizes. Since we expect the Berry phase to depend on the magnetization (that is not a sharp quantum number and can fluctuate), this is not surprising. Moreover, our results suffer from finite-size effects: the charge profile of the two vortices may overlap with each other and the charges accumulated at the edge of the finite system.  The results show good agreement with the theoretical prediction of φ π−flux − φ bare = π/2. For each system size, we have two different values of the number of particles N. We have q = 1/2 (q = −1/2 for the larger (smaller) of the two values of N.
As a benchmark for our method, we also performed calculations for the Berry phase accumulated by taking the vortex around a local flux π inserted in one triangular plaquette. The results are shown in Table. II. In this case, for a vortex of charge q, we expect a Berry phase of φ π−flux −φ bare = (2n+1)πq where n is an integer. With (without) an occupied vortex bound state the charge of the vortex was 1/2 (−1/2). Interestingly, for both cases we found a Berry phase close to +π/2, which is consistent with the above expression for n = 0 and n = −1 respectively (flux π is equivalent to flux −π). The benchmark above shows that while our numerical method is capable of reproducing established results, there is error of order a few percent in the finite-size numerical lattice calculation.
We now present direct measurements of n ↑ + n ↓ and n ↑ − n ↓ in a box of size 3.5 × 3.5 centered at the core of the moving vortex. The results are shown in Fig. 10. We see that the measured charge exhibits some fluctuations (a few percent) as a function of the position of the vortex core (on a single path and between different paths) due to the finite-size effects. The fluctuations of the magnetization are larger due to the nonconservation of spin. In particular when a vortex bound state is occupied, the local physics can strongly affect the magnetization of a vortex. Despite all these issues, the Berry phases we measured directly show good agreement with Eq. (16). In particular, if we average the fluctuating magnetization over the full cyclic path, the result is independent of the path. We can now compute the expected Berry phases from Eq. (16). The results are shown below and their agreement with Table I strongly supports the correctness of Eq. (16).

VIII. CONCLUSION
In summary, we have shown that integer quantum Hall systems, which emerge from the interplay of itinerant electrons and noncoplanar magnetic ordering, generically support topologically stable excitations with fractional charge and anyonic statistics. We showed through extensive numerical calculations that the statistical angle of these anyons has a simple relationship to their charge and magnetization.
The energetics of Z 2 vortices is similar to the usual Z vortices for S O(2) order parameter [39]: i.e., the energy of an isolated vortex scales logarithmically with the system size. Nevertheless, pairs of log-confined vortex pairs will appear due to thermal fluctuations at finite temperatures. In addition, inclusion of quantum spin dynamics may lead to an intriguing possibility that the quantum fluctuations transform the noncoplanar ordered state into a chiral spin liquid [48][49][50] at zero temperature. There, the fractionally charged Z 2 vortices discussed above may turn into deconfined point-like excitations.
Promising candidates materials which may exhibit this physics could be the systems of Na x CoO 2 type, which near x = 0.5 are known to have a noncollinear order, as well as anomalously large Hall response [51]. The fractional charge predicted in this work may be accessible through direct imaging of the local charge profile, as shown in Fig. 7, e.g., by scanning force microscopy. Also anyonic exchange statistics may have unusual consequences in real materials. Perhaps the most intriguing among them is the possibility of the Anyonic superconductivity [50,52,53]. When a system is doped away from the chiral Mott insulating state, it may energetically prefer to accommodate the carriers by creating vortices with intragap states. As we have just argued, such occupied vortex states are anyons, which, at finite density and low enough temperature, may go into a superconducting state [53]. Here we show the details of the derivation of Eq. (15) by explicitly writing out the terms appearing in two sides of this relationship. First, Eq. (14) for the three lattice vectors a i gives The Hamiltonian in the presence of the above δ 1 can then be written as which upon expanding to linear order in δ j gives H(k) = H 0 (k) + V(k) with H 0 (k) = −2 cos (k · a 2 ) e iπ/4 cos (k · a 1 ) + e −iπ/4 cos (k · a 3 ) e −iπ/4 cos (k · a 1 ) + e iπ/4 cos (k · a 3 ) − cos (k · a 2 ) , and The current operators can be obtained by differentiating H 0 (k) with respect to k, and are given by Inserting Eq. (A1) into Eq. (A4) and comparing with J x (k) and J y (k) above leads to Eq. (15):

Asymptotic flux pattern
As mentioned before, despite the relationship (15), the flux pattern is different than a localized flux π at the center of the vortex. Here, we explicitly show this intricate pattern. Far away from the vortex, the flux pattern has translation invari-type (a, b, c) δφ ance (to leading order) within local regions defined in Fig. 6 and can be easily calculated for a vortex with any axis of rotation l (so far we have only focused on l =ẑ). Consider a triangular plaquette with three magnetic moments S a,b,c, . If the effective flux through the plaquette is φ, we have the solid angle formula where | S a S b S c | indicates the determinant of a matrix A i j = S i ( j) (the jth component of S i ). It is easy to see that the numerator vanishes for any triangle in the unperturbed magnetic structure of Fig. 5(a), as it should for π/2 flux per triangular plaquette. Now if a moment S i is rotated by an angle θ i around l, we can show after some algebra that S a · S b → S a · S b cos δ ab + l · S a × S b sin δ ab + l · S a l · S b (1 − cos δ ab ). For small δ ab = θ a − θ b [see Eq. (14)], the leading correction comes from the second term. The leading correction to the flux then follows from the solid-angle formula [see Eq. (A7)], and the results are shown in Table III.
Given the intricate flux pattern, it is helpful to give a more microscopic derivation of the fractional charge using Eq. (15) and the Laughlin's argument. Basically, we want to compute the expectation value of the current flowing toward the vortex core (i.e., − J r ) in linear response. Instead of adiabatically inserting a local flux, however, in this case, we insert an intricate pattern of fluxes globally. In the region of Fig. 6, inserting this flux pattern corresponds to adiabatically turning on the perturbation V(k) of Eq. (15) from zero to its final value.
Let us consider a linear in time protocol with total time T as follows It is worth mentioning that it is common in the literature to add a term of the form e t H with → 0 + to the Hamiltonian in order to model adiabatically turning on a final perturbation H from t = −∞ to any finite t. Here, we use the above linear protocol, which is convenient for the case of quantum Hall response.
Note that, to invoke adiabaticity, we need the spectral gap to remain open during this process. We have checked numerically that the gap remains open if all the additional fluxes (due to the presence of the vortex) in each triangular plaquette is turned on linearly in time. The linear-response expression for the expectation value of an operator (in the present case J r ) at time t is given by where |0 is the initial state (in this case the ground state of H 0 in the absence of a vortex) and the "hat" notation represents an operator in the interaction picture with the bare Hamiltonian H 0 , e.g., To proceed, we diagonalize H 0 as follows: H 0 = k,a=1,2 ε a k γ † ak γ ak , where ε 1,2 k = ±2 cos 2 k · a 1 + cos 2 k · a 2 + cos 2 k · a 3 are the eigenvalues of H 0 (k) with corresponding orthonormal eigenvectors |1 k and |2 k . We can now write Eq. (A9) as where O mn (k) indicates m k |O(k)|n k . Summation over m, n, m , n gives 16 types of commutators. It turns out, however, that only the following two types of commutators have a nonvanishing contribution (i.e., reduce to γ † γ): and a similar commutator with 1 and 2 indices switched. Defining n m k = 0|γ † mk γ mk |0 , this leads to The integral over t can be simply done by integration by parts. We then need to integrate the resulting expression, which is independent of time, over time from 0 to T and over a circle of radius r for 0 < θ < 2π, which, putting back the factors of e and gives where the quantized Hall conductance σ xy is given by in the case of the integer Quantum Hall effect as in our system. All conductivities σ ab xy can be computed from the general expression (A12), but with the charge currents replaced by the appropriate charge or spin current J nm xa (k) and J nm yb (k). This expression is applicable to translationally invariant systems where momentum k is a good quantum number. For a unit cell of M sites, the momentum-space Hamiltonian H(k) can be generically written as a 2M×2M matrix. (The factor of two accounts for electron spin.) For each momentum k, we can diagonalize this 2M × 2M Hamiltonian and obtain its eigenvalues and eigenvectors: H(k) = ε m k |m k m k |, m = 1 . . . 2M. The eigenvalues ε m k give the 2M energy bands. If we have symmetries, as in the triangular lattice case discussed below, these bands can be degenerate. The matrix elements J nm xa (k) can be constructed explicitly using the eigenvectors |m k and |n k , where J(k) is an appropriate current operator written as a 2M × 2M matrix in the same basis as H(k).
The only ingredients for computing σ ab xy are then the 2M × 2M Hamiltonian H(k) and the corresponding 2M×2M charge and spin current operators J x,y a (k) for a = 0 . . . 3. With these ingredients, one can diagonalize H(k) to obtain the eigenvalues and eigenvectors, use the eigenvectors to construct the matrix elements of the current operators, perform the sum over m and n, and finally integrate the resulting expression over momenta k in the Brillouin zone.
Based on the above prescription, our main task is to write H(k) and J x,y a (k). We first choose an explicit tetrahedral magnetic ordering represented as in Fig. 5. As in the main text, we also assume an additional Zeeman field h in the z direction. We choose the following basis: , to write the Hamiltonian as H T = k Ψ † k H T (k)Ψ k , where the subscript T indicates the triangular lattice. We can then write where 3 (k) 2 (k) 1 (k) 0 2 (k) 3 (k) 3 (k) 2 (k) 0 1 (k) 2 (k) 3 (k) 1 with i (k) ≡ −2t cos(k · a i ). Charge (a = 0) and spin (a = 1, 2, 3) current operators can then be simply written in the same basis as An almost identical expression can be written for the kagomé lattice by choosing explicit components for the magnetic moments and a basis , c † 3↓k ) as shown in Fig. 11: (B5) The idea is to simulate the adiabatic motion of the vortex along a path by discretizing the path, and computing the many-body ground-state wave function of the system for the vortex core lying on each such discrete point along the path. Let us represent these wave functions by |Ψ i , with i = 1 . . . N + 1, and the periodic identification |Ψ N+1 = |Ψ 1 . We can then approximate the Berry phase Ψ|∂ s |Ψ ds by [12] which is a convenient expression for numerical calculations as it is explicitly gauge-invariant (each ground-state wave function obtained from exact diagonalization has an arbitrary U(1) FIG. 12. A typical path in a system of L = 40 with radius r = 0.2L. The discretization points are distributed at uniform angular coordinate, but the radial coordinate is shifted to maintain the minimal distance of 0.25 (in units of lattice spacing) from the lattice sites.
phase; however, that phase obviously drops out from this expression). The only underlying assumption for the validity of the above expression is 1 − | Ψ i |Ψ i+1 | 1, which requires having close discretization points so the overlaps | Ψ i |Ψ i+1 | are close to one.
For the case of hard-core vortices in a noncoplanar magnet, the wave functions can exhibit violent changes if the "moving" vortex core approaches very near a magnetic moment (i.e., a lattice site). Thus, in order for the Berry phase to con-verge faster, it is important to use more discretization points in regions of the path close to a lattice site, or, alternatively distort the path in the vicinity of sites so as to maintain a minimum distance from them. Here, we choose the latter approach with a typical path shown in Fig. 12. We first considered equally spaced discretization points (core of the moving vortex) on a circle centered at the other (fixed) vortex. For each discretization point we then find the distance to the closest lattice site and, if necessary, distort the path by shifting the point along the radial direction to the distance of 1/4 away from the site. This simple trick results in large | Ψ i |Ψ i+1 | overlaps.
To obtain reliable Berry phases from such numerics, we must have convergence in the number of discretization points (which requires large overlaps of consecutive wave functions). Also, due the presence of edge modes for open boundary conditions, we must be careful about level crossings at the Fermi level. We work at a chemical potential inside the gap but as the energy of the intragap states depends on microscopic details (like the current location of the "moving" vortex core inside a unit cell), for some trajectories at fixed particle number, the last filled level can switch from an edge mode to a bound state, rendering the results unreliable. This can be easily checked a posteriori by keeping track of the energies of the intragap levels as the vortex moves along the path.