Non-Abelian statistics as a Berry phase in exactly solvable models

We demonstrate how to directly study non-Abelian statistics for a wide class of exactly solvable many-body quantum systems. By employing exact eigenstates to simulate the adiabatic transport of a model's quasiparticles, the resulting Berry phase provides a direct demonstration of their non-Abelian statistics. We apply this technique to Kitaev's honeycomb lattice model and explicitly demonstrate the existence of non-Abelian Ising anyons confirming the previous conjectures. Finally, we present the manipulations needed to transport and detect the statistics of these quasiparticles in the laboratory. Various physically realistic system sizes are considered and exact predictions for such experiments are provided.


Introduction
A striking feature of topological phases of matter is that they can support anyons. These are quasiparticles with statistics different from bosons or fermions. The statistical behavior of anyons is demonstrated by their adiabatic exchange, which causes non-trivial evolution in their quantum state. For Abelian anyons the evolution is given by a phase factor. The presence of non-Abelian anyons gives rise to degeneracy in the energy spectrum and the evolution is described by a unitary matrix acting on the degenerate states. In general, anyons are known to exist in different varieties distinguished by their characteristic statistics [1]. Therefore, the explicit demonstration of the statistical behavior is essential for the unique characterization of a topological phase. Such phases are of great interest due to the possibility of realizing anyons in a physical system and due to their potential for technological applications. In particular, topological quantum computation employs the anyonic statistics for performing error-free quantum information processing [2].
The best known many-body system conjectured to support non-Abelian statistics is the fractional quantum Hall liquid [3,4,5]. Other proposals include the p-wave superconductor [5,6] as well as various lattice models [7,8,9]. These systems are either tailored to identically support non-Abelian statistics and have complex physical realizations, or they can be described by simple Hamiltonians, but their statistical behavior is based on indirect arguments. In particular, for the fractional quantum Hall states they rely on properties of trial wave functions [10,11,12], whereas for the lattice models explicit calculations have not been previously attempted. Although the indirect arguments are sound, direct calculations of the statistics are crucial to resolve any ambiguities, to address physical realizable finite-size systems and to provide exact predictions for the experiments.
Here we demonstrate how to directly calculate the non-Abelian statistics for a class of exactly solvable models. By applying the Berry phase technique [10] to the Kitaev's honeycomb spin lattice model [9], we calculate the evolution associated with an adiabatic exchange of quasiparticles. This is performed using exact eigenstates of a 360 spin system. We obtain a unitary matrix that corresponds to the statistics of the conjectured non-Abelian Ising anyons. Together with the fusion rules of these anyons [13,14], this conclusively demonstrates the non-Abelian character of Kitaev's model, thereby confirming the conjectured behavior. Further, we present a scheme for creating, transporting and characterizing the anyons that could be used in the proposed physical implementations [15] and provide exact predictions for a physically realistic range of the model's parameters.
is topologically equivalent to two unlinked loops.

The honeycomb lattice model
Kitaev's model [9] comprises of spin-1/2 particles residing on the vertices of a honeycomb lattice. The spins interact according to the Hamiltonian where J ν ij are positive nearest neighbor couplings on links (i, j) of type ν (see Figure 1(a) for link labeling). The second term is an effective magnetic field with positive next-tonearest neighbor couplings K ijk , such that every plaquette p contributes the six terms The enumeration of the sites is shown in Figure 1(a). The Hamiltonian has the symmetry [H,ŵ p ] = 0, whereŵ p = σ x 1 σ y 2 σ z 3 σ x 4 σ y 5 σ z 6 are plaquette operators whose eigenvalues w p = −1 are interpreted as having a vortex on plaquette p. We represent the spin operators as [9,13]. Subsequently, the Hamiltonian takes the form Here J ij andû ij are shorthand notations for J ν ij andû ν ij when (i, j) is an ν-link. Since the mapping to Majorana fermions doubles the size of the Hilbert space, the eigenstates of the original Hamiltonian (1) are subject to the constraint which follows from the operator identity σ Since [H,û ij ] = 0, the Hilbert space splits into sectors each labelled by u, a certain pattern of eigenvalues u ij . The configurations u can be understood as a classical Z 2 gauge field with local gauge transformation operators D i . Consequently, the plaquette operatorsŵ p = (i,j)∈pû ij can be identified with gauge invariant Wilson loop operators, whose patterns of eigenvalues label the physical sectors of the model. Fixing the gauge field configuration u gives then a particular vortex configuration. Throughout this paper we use J ν , K and u without indices to denote global configurations of these local quantities and use indices only when referring to their local values. For instance, K = a means that K ijk = a for all i, j and k.
Consider the system of 2MN spins on a torus and assume u to be fixed such that it creates the four vortex configuration shown in Figure 1 and ǫ k are the corresponding positive eigenvalues. In [13] it was shown that when the system is in the non-Abelian phase (J x = J y = J z = 1, K > 0), the presence of 2n well separated vortices gives rise to n zero modes (ǫ k ≈ 0, for k = 1, . . . , n). Importantly, these are separated from the rest of the fermionic spectrum by a finite energy gap. In our case of four vortices this implies fourfold ground state degeneracy arising from a pair of zero modes that can be either occupied or empty where α 1 , α 2 = 0, 1 and | gs = M N k=1 b k | φ is the ground state. For convenience we choose the reference state such that b † k | φ = 0. Numerical diagonalization of A, (see (2)), gives 2MN eigenvectors ψ ± k satisfying the double spectrum Aψ ± k = ±ǫ k ψ ± k , where ǫ k coincide with the positive eigenvalues of the diagonalized Hamiltonian. We construct a representation of the two degenerate ground states | Ψ 10 and | Ψ 01 as where α = 1, 2, respectively, and ε k,...,l is the fully anti-symmetric tensor of rank MN −1.
In general, such states are too large to be stored in a computer, because their number of elements grows exponentially with the system size. However, the inner product of two such vectors, each depending possibly on some parameters t and t ′ , can be efficiently calculated and is given by where

The Ising anyon model
It has been conjectured that Kitaev's model supports the Ising anyon model [9,3,5]. This model has three types of particles: 1 (vacuum), ψ (fermion) and σ (non-Abelian anyon). In [13] these are identified with the ground state, the fermion modes b † and the vortices, respectively. The non-trivial fusion rules are given by ψ × ψ = 1, ψ × σ = σ and σ × σ = 1 + ψ. The last fusion rule implies that there is a degree of freedom associated with the different ways a number of σ's can fuse when their total anyonic charge is fixed. Taking the four σ particles to fuse to a ψ, this fusion degree of freedom is encoded in the two dimensional fusion space, V ψ σ 4 . Its basis can be chosen to be the states associated with the two distinct pair-wise fusion channels: In [13] the number of intermediate ψ's is identified with the number of occupied zero modes. Hence, a suitable basis is given by the states {| Ψ 1 , | Ψ 2 } (see (5)). The braid operator, R, describes the statistics of the σ anyons. In particular, the monodromy operator, R 2 , corresponds to one particle encircling another clockwise. On the basis (7) the monodromy of two σ's that belong to different pairs is given by

Non-Abelian statistics as a holonomy
When z 1 and z 2 are the coordinates of the σ anyons, their statistics is given by the transformation of the wave function under their permutation, i.e. ψ(z 1 , z 2 ) = Uψ(z 2 , z 1 ) with U being the characteristic statistical phase or matrix. In real physical systems the permutation of the coordinates corresponds to adiabatically transporting the anyons such that their positions are swapped. When the positions are swapped twice, i.e. a particle winds around the other along a suitable chosen path, the statistics corresponds to the accumulated wave function evolution, which is given by the Berry phase, or the holonomy [10,16]. For bosons and fermions this is always trivial, with non-trivial evolution being a sign of anyonic statistics. We demonstrate the statistics of σ anyons by adiabatically transporting a vortex around another. Consider a Hamiltonian H(λ) with n-fold degeneracy {| Ψ α (λ) |α = 1, . . . , n} that depends on some parameters λ. When we adiabatically vary λ along a closed path C, the evolution of the degenerate subspace is given by the holonomy and P denotes path ordering in λ. To simulate the vortex transport, we discretize the path C into T infinitesimal intervals of length δλ with λ(t) denoting the control parameter value at step t. It follows that the holonomy takes the form i.e. in the limit δλ → 0 it is given by the ordered product of projectors onto the ground state space at each step t.
We evaluate the evolution in the fusion space V ψ σ 4 . The basis states (5) are not symmetrized under gauge transformations (3). Nevertheless, their holonomy coincides with the holononomy of symmetrized states when C is a loop in both the space of four vortex and gauge field configurations. This is due to the orthogonality of states belonging to different sectors of u. A suitable path is illustrated in Figure 1(a), where the path C (dashed lines) is split into four parts. Different ordering of these parts corresponds to the topologically inequivalent paths C l and C o given in Figure 1(b) and (c), respectively. Neither path spans any area and hence all contribution to the holonomy is topological. Since u is a static background field, we need to introduce classical control parameters to physically implement the transport. Assuming local control of J ij and K ijk on all links, we see from (2) that the simultaneous sign change of these quantities on link (i, j) is equivalent to changing u ij → −u ij . This either generates a vortex pair or transports a vortex through the link (i, j). In our simulation this is performed in S infinitesimal steps. Taking λ = (J, K) and assuming T to be a sufficiently large, the discrete holonomy (9) for the degenerate states (5) is well approximated by where we have used the inner product (6) ‡. Therefore, the holonomy can be evaluated by diagonalizing the Hamiltonian at each step t and multiplying together the inner products of the eigenstates from successive steps according to (10). We perform this for the three parametrizations shown in Table 1. Table 1. Three parametrizations (i), (ii) and (iii) for which the holonomy is evaluated.
Here T = 8S(d + 1) and the number of spins is 2M N = 8(d + 2)(2d + 3). S has been increased in (iii) to suppress accumulation of discretization errors due to longer path.
Since the spectrum varies slightly with t during the braiding process, we define the minimal fermion gap, ∆, and the maximum energy splitting between the two ground states, δ, by respectively, where ǫ t k is the kth eigenvalue at step t. These are plotted in Figure 2, where we observe that both the fermion gap and the level of degeneracy improve as K and d increase. Under the adiabatic approximation the holonomy corresponds to the exact time evolution when ∆ ≫ δ and δ → 0. To physically accommodate these ‡ The freedom in changing the basis at each step t of the path C gives rise to an accumulated unitary matrix, M , making the result of the adiabatic evolution to be in general given by M Γ C [12]. Here we choose | Ψ α (t = 0) = | Ψ α (t = T ) which gives M = 1 1.  Figure 2. The minimal fermion gap ∆ (--) and the maximum energy splitting between the ground states δ, (11) (----) as functions of K for parametrizations (i) (• ), (ii) (⊓ ⊔) and (iii) (♦) given in Table 1. The fermion gap grows linearly and the degeneracy improves with increasing K for all parametrizations. The fermion gap is relatively insensitive to the vortex separation, whereas the degeneracy improves when the vortices are further apart.
conditions in a finite size system, the vortex transport should be fast enough compared to δ for the states | Ψ α , α = 1, 2 to appear as degenerate, but slow enough compared to ∆ so that no fermionic excitation is produced. We see from Figure 2 that δ ∆ 10 −2 for (iii) when K 0.07. This region can support the adiabaticity conditions and hence we take K ≈ 0.07 as a lower bound for identifying a stable topological phase.
To quantitatively study the holonomy, we introduce a fidelity measure for a target matrix U and a test matrix V as For U, V unitary 2 × 2 matrices we have that s(U, V ) = 1 if and only if U = V , while in general s(U, V ) ≤ 1. Let Γ C l be the numerically obtained holonomy with off-diagonal elements re iθ , 0 ≤ r ≤ 1. After fixing the gauge §, we evaluate the unitarity measure, s(1 1, Γ C l Γ † C l ), Figure 3(a), and the two different fidelity measures of the holonomy: s(|R 2 |, |Γ C l |) = r (measure of off-diagonality that characterizes R 2 ) ands(R 2 , Γ C l ) = 1 2 [s(R 2 , Γ C l ) + 1] = 1 2 [r cos( π 4 + θ) + 1] (the total fidelity), Figure 3(b-d). Here |U| denotes a matrix U with its elements replaced by their absolute values.
First, we observe that the unitarity measure is above 98% for all parametrizations when K 0.10, which we take as an upper bound for identifying a stable topological phase. For (i) we obtain no significant off-diagonality due to the small size of the system. However, for (ii) the holonomy is predominantly off-diagonal (e.g. r > 0.9) for 0.02 K 0.04, and for (iii) for 0.02 K 0.09. The total fidelity,s, accounts also for the overall phase and can distinguish between the Ising (s = 1) and SU (2) Table 1. The measure of off-diagonality, s(|R 2 |, |Γ C l |) (--), and the total fidelity,s(R 2 , Γ C l ) (----), as a function of K for the parametrizations (b) (i) (• ), (c) (ii) (⊓ ⊔) and (d) (iii) (♦). Based on unitarity and the energy gap behavior, we expect a stable phase in the area 0.07 K 0.10 bounded by the dashed vertical lines.
anyon models whose monodromies only differ by an overall phase factor, e −iπ/2 [1]. We observe that for (ii) there is a small region around K ≈ 0.02 and for (iii) there is a wider region, 0.08 K 0.10, wheres > 0.9. The maximum fidelities are given by 0.981 and 0.991, respectively. Parametrization (iii) also has a region 0.02 K 0.05 wherē s ≈ 1 2 with error ±10 −1 . However, we disregard this regime, because such a region does not exist for the smaller system (ii) and it lies outside the domain which we consider as a stable topological phase. Further, we check for all parametrizations and all K that Γ Co ≈ 1 1 with error less than 10 −2 , that for K = 0 the holonomy vanishes and that Γ C −1 l = Γ † C l when the direction of braiding is reversed.

Braiding and detection of the non-Abelian statistics in a laboratory
Since our calculation involves only the experimentally accessible parameters J and K, it translates directly to how one could physically implement the creation and transport of anyons in the laboratory. In particular, in the optical lattice proposal of Micheli et al. [15], the vortex transport would correspond, given sufficient site addressability, to the local adiabatic inversion of magnetic field as well as of the couplings J through the introduction of suitably tuned lasers. Also, as the energy of the system is known to depend on the zero mode populations [13], the effect of braiding can be detected through spectral means. As the monodromy swaps these populations between the vortex pairs, the energy behavior of the system will be different when the vortices from a single pair are brought close together before and after the braiding. Detecting this energy shift reveals the non-Abelian statistics. By evaluating the holonomy, we were able to identify a range of the model's parameters where the simulation approximates well the exact time evolution of a physical system and where the statistics corresponds to the Ising anyons. As expected, larger systems exhibit the predicted statistics with higher fidelity. The required magnitude of K, however, is larger than anticipated (K ≈ 0.1 for parametrization (iii)). In the original work [9], the three-body term appears in third order perturbation theory when one considers a general Zeeman term (h ν i σ ν i ) as a perturbation. In our normalization the expansion is valid when h 2 ≪ 1. Since K ∼ h 3 , for K = 0.1 one can estimate h 2 ≈ 0.2, which clearly does not satisfy the criteria. Therefore, in order to introduce the three-body terms into the Hamiltonian perturbatively, such as by adding a small magnetic field in the optical lattice proposal [15], one needs to consider larger systems. On the other hand, were the three-body terms engineered [17], our calculation provides exact predictions for braiding experiments in such systems.

Conclusions
In summary, we formulated a method to directly study non-Abelian statistics in exactly solvable lattice models whose ground state admits representation as a Slater determinant. By applying it to Kitaev's model, we identified finite regions of the couplings, where the non-Abelian statistics corresponds to Ising anyons. This confirms the previous conjectures for the presence of a non-Abelian topological phase. Finally, we proposed a scheme for the implementation and detection of non-Abelian statistics in the laboratory. Such an experiment would be an important step towards the physical realization of topological quantum computation. It is an interesting topic for future research to study whether the holonomy can be used as an order parameter for the topological phase when the system is subject to perturbations.