Rigorous Calculations of Non-Abelian Statistics in the Kitaev Honeycomb Model

We develop a rigorous and highly accurate technique for calculation of the Berry phase in systems with a quadratic Hamiltonian within the context of the Kitaev honeycomb lattice model. The method is based on the recently found solution of the model which uses the Jordan-Wigner-type fermionization in an exact effective spin-hardcore boson representation. We specifically simulate the braiding of two non-Abelian vortices (anyons) in a four vortex system characterized by a two-fold degenerate ground state. The result of the braiding is the non-Abelian Berry matrix which is in excellent agreement with the predictions of the effective field theory. The most precise results of our simulation are characterized by an error on the order of $10^{-5}$ or lower. We observe exponential decay of the error with the distance between vortices, studied in the range from one to nine plaquettes. We also study its correlation with the involved energy gaps and provide preliminary analysis of the relevant adiabaticity conditions. The work allows to investigate the Berry phase in other lattice models including the Yao-Kivelson model and particularly the square-octagon model. It also opens the possibility of studying the Berry phase under non-adiabatic and other effects which may constitute important sources of errors in topological quantum computation.


Introduction
Quantum statistics is intimately related to the adiabatic exchange of indistinguishable particles. Exchanging two particles twice results in a loop trajectory which in three dimensional space can be smoothly contracted to a point, equivalent to no trajectory. The particles' wavefunction thus remains unchanged after two subsequent exchanges, and after one exchange can transform either in a symmetric or an antisymmetric fashion, giving Bose-Einstein and Fermi-Dirac statistics respectively. In two dimensional space such a contraction of the particles' trajectory is impossible. This gives rise to a different type of quantum statistics. Algebraically, adiabatic exchange operators correspond to the elements of the permutation group S N for three and higher dimensional systems, and to the elements of the braid group B N for two dimensional systems [1,2]. Both groups are formed by N − 1 generators τ 1 , ..., τ N −1 , obeying the constraints where the generator τ i interchanges the two particles at the positions i and i + 1. The quantum statistics then arises from the unitary irreducible representations (irreps) of these groups [1]. The group S N possesses two one-dimensional irreps which correspond to bosonic and fermionic statistics. Its higher dimensional irreps can be replaced by bosonic and fermionic statistics when a hidden degree of freedom is introduced [3].
On the other hand, one-dimensional representations of the braid group B N can be any phase e iθ where θ ∈ [0, 2π) (hence the name anyonic) [4,5]. More interestingly, the braid group also permits multi-dimensional unitary irreducible representations which give rise to non-Abelian statistics. Any exchange of particles then leads to a unitary rotation of a state vector of the system within a D-fold degenerate ground state. The degree of degeneracy D depends only on the presence of N well-separated identical particles. As this is solely linked to the topology of the underlying configuration space, braiding the particles is the only way to induce nontrivial operations within this ground state subspace. Consequently, the system is immune to any local perturbations or fluctuations as long as these do not exceed the spectral gap which separates the rest of the system's spectrum from this decoherence free subspace [6]. This capability to implement unitary operations within an intrinsically fault-tolerant framework offers promising applications in quantum information processing, specifically topological quantum computation [6,7].
A larger body of research results suggests that non-Abelian anyons are physically realized as localized quasiparticle excitations of many-body systems. These for example include the systems that manifest the fractional quantum Hall (FQH) effect [8,9], p x + ip y superconductors [10] and the Kitaev honeycomb spin lattice model [11] (with its proposed realizations [12,13]). Theoretical studies show that all these systems can be effectively described by a similar topological quantum field theory (the Ising and the related SU(2) 2 Chern-Simons theory) which can be characterized by three particle types labeled as 1 (i.e. vacuum or trivial topological charge), ǫ (fermion) and σ (the non-Abelian anyon). These satisfy certain fusion and braiding rules which will be specified later. The experimental study of non-Abelian anyons is indeed of fundamental importance and so far the experimental observations have yielded encouraging results, though it is to be said that final verification of non-Abelian statistics remains a great challenge [14].
Physically, anyonic statistics arise from the evolution of the system under adiabatic interchanges of these quasiparticles. According to the adiabatic approximation, a physical system remains in its instantaneous eigenspace if a given perturbation is acting on it slowly enough and if there is a gap between the corresponding eigenvalue and the rest of the Hamiltonian's spectrum. When these perturbations draw a smooth and closed trajectory C(λ) in the parameter space, the unitary evolution of the system in the n dimensional eigenspace is given by the Berry phase (matrix) B(C) B(C) := P exp i A(λ)dλ (4) where P denotes the path-ordered integral and where |Φ l (λ) and |Φ k (λ) are the eigenstates of the system's Hamiltonian at the value of the parameter λ.
In this paper, we directly evaluate the non-Abelian statistics of the Ising anyons of the Kitaev honeycomb model. In particular, we numerically calculate the non-Abelian Berry phase (matrix) which governs the evolution of the system under the adiabatic exchange of two σ-particles (vortices) of the Kitaev honeycomb model. This work can be seen as an accurate test of the non-Abelian statistics in the Kitaev model which offers applications in the context of topological quantum information processing and computation. Moreover, it provides a direct way to study the non-Abelian statistics in lattice models with a quadratic fermionic Hamiltonian and as such it complements similar efforts carried out in the context of continuous systems [15][16][17].
The non-Abelian Berry phase calculations in the Kitaev model have been a subject of study by Lahtinen and Pachos [18]. They developed an interesting technique for inducing the vortex motion in the Kitaev honeycomb lattice model which we have utilized in the present work, though within a different solution of the model. While the previous studies established the non-Abelian nature of the statistics, the results on both the exact form of the braid matrix and on the exponential convergence with vortex separation were not conclusive.
We establish these results rigorously for much larger vortex separations and to a very high accuracy, but we would like to emphasize that our approach goes beyond a mere technical improvement. Our calculations rely on the solution of the Kitaev model which was presented recently by Kells et al. [19]. This solution employs the Jordan-Wignertype fermionization in the exact effective spin-hardcore boson representation of the model and uses no redundant degrees of freedom, thus allowing us to work directly with physical eigenstates of the system. This allows us to calculate the Berry phase associated with braiding vortices at the minimal distances for up to nine plaquettes. The simulation also reaches a very high degree of accuracy as measured by the Frobenius distance between the Berry matrix obtained from our simulation and the exact Berry phase from the effective field theory. The accuracy of our calculations increases exponentially fast with the vortex distance, achieving results which are characterized by errors on the order of 10 −5 and lower.
We thus present in this work a very accurate technique for the non-Abelian Berry phase calculation. The accuracy of the calculations allows us to see the exact dependence of the simulated Berry phase on the details of the model, like the splitting of the ground state level which is intrinsic to any finite system. It also shows the dependence on an exact implementation of the braiding operations, potentially allowing for an analysis of non-adiabatic effects, similar to that provided quite recently in a different context in [20]. The importance of these effects derives from the fact that they will likely be a crucial source of errors for any topological quantum information processing and computation. Possible applications thus extend to modeling and implementation of quantum information protocols whose reliability can be tested under various effects; for example, disorder. Naturally, the first step in this potentially fruitful story is the demonstration of highly accurate and sufficiently large scale direct simulation of the Berry phase in the Kitaev model, as presented in this manuscript.
The results we present show strong agreement with the statistical properties of Ising anyons derived from the effective theory. One can naturally wonder what the meaning of calculating the Berry phase is if we already know it from the relevant effective field theory. The point here is that the effective theory gives us nearly no ground to test the stability of the Berry phase under the effects mentioned above or to make any conclusions about its implementation under (more) realistic conditions and about what pitfalls we should expect in such situations. Moreover, the accurate calculations like those presented here provide important predictive power in the analysis of topological phases in other lattice models including, for example, the Yao-Kivelson model [21] and the square-octagon model [22] which exhibits a kaleidoscope of topological phases including the Ising and SU(2) 2 phases. We believe that this is highly relevant to implementation of quantum information processing in Majorana fermion systems [23].
The paper is organized as follows. In Section 2, we review the exact solution of the Kitaev model on a honeycomb lattice by explicitly describing its eigenstates based on the method presented in [19]. We first define Jordan-Wigner types fermions in the honeycomb lattice and then represent the original Hamiltonian using these fermions. The resulting Hamiltonian is in a quadratic fermionic form which can be solved exactly. Then in Section 3, we discuss how to smoothly move the vortices within our solution of the model and investigate the adiabaticity of the anyonic motion. When the anyons follow a cyclic path adiabatically, the evolution of the system is governed by the Berry matrix; the numerical method for calculating this matrix will be presented in Section 4. The presented method can be applied to the Berry matrix of any system having a quadratic fermionic Hamiltonian. In Section 5, we discuss the results of the numerical calculation by comparing them with the expected statistics and present an error analysis of the numerical calculations. Vertex coloring emphasizes the bipartite lattice structure. The operators W p , S q,• , L x and L y are defined as products of single-and two-body terms. The Kitaev honeycomb model is a spin-1/2 lattice model in which spins are located on the vertices of a honeycomb lattice (see Figure 1(a)), and it has the following Hamiltonian:

The Kitaev Honeycomb Model
where i and j are the position indices of the spins, J α i,j , α = x, y, z, are the coupling coefficients of the two-body interaction operator K α i,j = σ α i σ α j on the link (i, j), and the σ α i are the Pauli operators. The model is exactly solvable and contains three equivalent gapped A phases for parameters satisfying J x > J y + J z , J y > J x + J z or J z > J x + J y , and a gapless B phase for the other values of the parameters. Furthermore, adding to the Hamiltonian (5) a term which breaks the time-reversal and parity invariance of the model opens a spectral gap in the B phase and allows the realization of non-Abelian anyons of the Ising type. This additional term, defined as where p indexes honeycomb plaquettes, represents the time-reversal and parity invariance breaking effect of a weak magnetic field [11] as it emerges from perturbation theory on the third order. The coefficients κ l p of the effective term l at the plaquette p (Figure 1(a)) are related to the magnetic field as κ ∼ hxhyhz In what follows we will consider only the effective magnetic field (6).
The model has a commuting set of plaquette operators for each hexagon p which also commute with the total Hamiltonian H tot = H + V. The eigenvalues of W p correspond to whether the plaquette p is occupied by a vortex (−1) or not (+1). The vortices carry unpaired Majorana modes for odd values of the Chern number ν [11]. For translationally invariant configurations of the Kitaev honeycomb model, it is found that ν = ±1 depending on the direction of the magnetic field [11]. Majorana modes exhibit non-Abelian statistics [24] corresponding to σ-particles of the Ising anyons. We will discuss their properties in more detail in Section 5. The model can be solved by various fermionization techniques, but here we will use the solution introduced in the paper [19] for our purposes. This solution has the advantage of giving the eigenstates of the system explicitly whereas it is not practically possible to do that in Kitaev's original solution. The original solution maps the spin degrees of freedom of the model to Majorana fermions. This requires each spin degree of freedom to be embedded in an extended Hilbert space of four dimensions. Obtaining physical states then requires projections from the eigenstates of the extended Hamiltonian which is hard to achieve in practice and thus limits the extent of numerical calculations. For the sake of self-completeness and clarity for further arguments, we start with a brief discussion of the solution we will rely on.

The exact solution of the model
Here we focus on N x ×N y lattices on a torus where N x and N y are the numbers of z-links in the n x and n y directions respectively, as in Figure 1(a), where the dotted lines define a 4-by-4 lattice whose opposite sides are identified. Let us label the z-links by q = (q x , q y ) with respect to the z-link at the origin O = (1, 1). Spins are denoted by either empty • or full • circles, reflecting the bipartite structure of the underlying honeycomb lattice.
Periodicity on the torus imposes p W p = 1 on the plaquette operators and gives rise to homologically nontrivial loop symmetry operators L x and L y (see Figure 1(b), (c)) which have ±1 eigenvalues and commute with all W p and the total Hamiltonian [25]. Therefore, for a 2N-spin system on a torus, there are N + 1 independent operators which split the total Hilbert space of the system into 2 N +1 different 2 N −1 -dimensional subspaces.
Before we define the fermions on the lattice, let us define the string operators between an arbitrary location on the lattice q = (q x , q y ) and the origin O = (1, 1) (see Figure 1(a)) as S q,• := S y S x σ x q,• and S q, • to the z-links and x-links of the interval [O, (q x , 1)) respectively and similarly S y is the successive applications of σ z • σ z • and σ y • σ y • to the z-links and y-links of the interval [(q x , 1), (q x , q y )) respectively. Note that S x = I when q x = 1 and S y = I when q y = 1. S q,• and S q,• commute with L y and all plaquette operators except the one located on the left of the origin (i.e. W (Nx,1) ) and L x , with which they anti-commute.
By using string operators, fermionic creation and annihilation operators are defined on each z-link q as It is not difficult to show that they satisfy the fermionic anti-commutation relations Because the c-fermions and the string operators have the same commutation and anticommutation relations with the plaquette and the loop operators, the quadratic forms of fermionic operators ( e.g. c q c † q ′ for any q and q ′ ) commute with L x , L y and W q for all q. In other words, for a 2N-spin system on a torus, the even fermionic states span the 2 N −1 dimensional subspaces of each {W p , L x , L y } configuration. Therefore, only states having an even number of c-fermions are realized.
By using these fermions, we write H in Equation (5) as where X and Y are defined (see Figure 2) as The transformation of a honeycomb lattice into a square lattice by contracting the z-links to a point. (b) and (c) show X (2,4) , Y (3,4) respectively for a 4 × 4 square lattice whose opposite sites are identified.
Similarly, the fermionic representation of the P l q terms of V (6) reads Notice that the total Hamiltonian H tot = H + V of the Kitaev model is quadratic fermionic. The quadratic fermionic Hamiltonian H of any system with N fermions must be of the form where ξ is Hermitian and ∆ is antisymmetric, and can be rewritten in the following way [26] where is a Hermitian matrix that can be diagonalized as whose columns correspond to eigenvectors of M. The matrix M D : E is a diagonal matrix with positive entries. These are placed in an increasing order E 1 < ... < E N as a convention [26]. By replacing M in Equation (12) Since T is a unitary matrix, it is possible to define new sets of fermions These two definitions are compatible, and therefore, the Hamiltonian H can be rewritten in a free fermionic form as follows: For the total Hamiltonian H tot of the Kitaev model, M tot -the analog of the matrix (13) -is given in terms of X and Y ; therefore it is diagonalized separately for each X and Y configuration (i.e. {W p , L x , L y } configuration).
On the other hand, the eigenstates of the quadratic fermionic Hamiltonians can be explicitly written by using the Bloch-Messiah theorem [27,28]. According to this theorem, any unitary matrix of the form of T in Equation (14) can be decomposed into three matrices of very special form: where D and C are unitary matrices and U and V are real matrices of the general block-diagonal form where Z and I are the zero and identity matrices of the same size respectively, and where u i and v i are positive real numbers. By using Equation (17) in Equation (15), we have Here, D is used to define new operators a † and a: Then there is a special Bogoliubov transformation This distinguishes the "paired" levels (u p > 0; v p > 0) where (p, p) are defined by U i and V i (18)) from the "blocked" levels (v m = 0; u m = 1) which are either occupied (v i = 1; u i = 0) or empty. Finally, a linear transformation of the α † and α by the unitary matrix C gives: For a general quadratic fermionic Hamiltonian, the ground state wavefunction is defined as a non-zero wavefunction |φ such that β k |φ = 0 for all k. It can be easily verified that the following wavefunction satisfies these criteria: where |− represents the vacuum of c-fermions. In the honeycomb model, the vacuum |− belongs to the X,Y -configurations for which M tot is diagonalized. However, because only the even fermionic configurations are allowed for each X,Y -configuration, when the number of elements in the first product (i.e. i-part) of (19) is odd, the ground state |Φ of the model is the next excited state |Φ = β † 1 |φ where β † 1 is the minimum energy fermion. In that respect, the number of a † i determines the fermionic parity of the system. For future reference, it is important to associate every eigenstate of the system with a T = U V * V U * matrix such that the Bloch-Messiah representation Equation (19) represents that eigenstate. For an excited state β † i |φ , this can be done by exchanging the roles of β i (i) and β † i (i) which manifests itself as the exchange of the ith and (N +i)th columns of T ; we denote this matrix by T ′ . Thus the ground state |φ ′ of T ′ is equal to β † i |φ . A straightforward generalization of this method to the other excited states associates any eigenstate with a particular T matrix.

Adiabatic Motion of Vortices and the non-Abelian Berry Phase
In this Section, we discuss how to move vortices from one plaquette to another, then describe the evolution of the system under an exchange of vortices.
By changing the sign of the relevant coupling coefficients J and κ, we can emulate the fermionic spectrum relevant to any X,Y -configuration starting from the trivial X, Yconfiguration (i.e. X q = 1 and Y q = 1 for all q) which we call the reference X, Yconfiguration.
This approach, first introduced in [18], allows us to consider the Hilbert spaces of different X, Y -configurations connected by the coupling coefficients J and κ.
In this way, we can also simulate the creation, annihilation and motion of the vortices. For example, consider the {W p , L x , L y } configuration (see Figure 3(a)) L x = L y = −1 and W q = 1 for all q except W q ′ = W q ′ − ny = −1; where q ′ = (q ′ x , q ′ y ) such that q ′ y = 1. This configuration gives Y q = 1 and X q = 1 for all q except X q ′ = −1. The fermionic spectrum of this configuration can be achieved from the reference X, Y -configuration by using the negative values of the set of coefficients (10) and (11)). In other words, changing the sign of can be considered as the creation of two vortices from the vacuum (i.e. the reference X, Y -configuration with positive J and κ values).
For an analogous configuration shown in Figure 3(b), vortices on the plaquettes q ′ and q ′ − n x for q ′ y = N y and q ′ x = 1 can be created from the reference X, Y -configuration by changing the sign of the set of coefficients [Jκ] y q ′ := {J y q ′ , κ 5 q ′ − nx , κ 6 q ′ , κ y q ′ }. Generally, changing the sign of [Jκ] x q alters the vorticity of two plaquettes W (qx,qy) and W (qx,qy−1) when q y = 1, and may be seen as creation, annihilation or motion of vortices, depending on the initial vorticity of the plaquettes. A similar effect can be achieved on the plaquettes W (qx,qy) and W (qx−1,qy) for q satisfying q y = N y and q x = 1 when we change the sign of [Jκ] y q . We point out here for completeness that it is also possible to simulate the vortices by changing the [Jκ] x q when q y = 1 or by changing [Jκ] y q for q satisfying q y = N y or q x = 1 by carrying the sign to some of the c-fermions. However we will use the previous approach to exchange vortices as it is sufficient. From now on we will use the term vortex only for these simulated vortices. In Figure  4 a template configuration is shown with four vortices which are created from the vacuum by changing the sign of [Jκ] x of the x-links between A and B, and between C and D (colored by red in Figure 4). In this paper, we will work with several configurations similar to Figure 4 with different sizes as: N x = 3d + 1, N y = 2d when d is odd and N x = 3d, N y = 2d when d is even for d = 1, ..., 9, where d is the minimum distance between the vortices. Note that all configurations are even-by-even; the other configurations (odd-by-odd, even-by-odd etc.) will be studied in future. For all these configurations, the coupling coefficients for the vacuum are identical for all plaquettes: J x q = J y q = J z q = J = 1 for all q. The strength of the effective magnetic field κ for the vacuum is also the same for all the plaquettes q, κ x q = κ y q = κ 5 q = κ 6 q = κ, and is taken from two sets which differ in magnitude: (i) κ = l × 0.01 where l = 1, ..., 5, and (ii) κ = l × 0.05 where l = 2, ..., 8.
To swap the position of the vortices B and C -see Figure 4 -we need to adiabatically move the vortices along the paths indicated by the arrows 1, 2, and 3. This requires us to slowly change the sign [Jκ] of the links which are intersected by the paths. For the large values of κ, i.e. the values from the set (ii) above, all these configurations have a (nearly) two-fold degenerate ground state which is separated by a gap from the rest of the Hamiltonian's spectrum. This can be seen from Figure 5(a) which shows the average of the ratio of the splitting of the ground state degeneracy and of the gap to the first excited state G 2,1 (λ)/G 3,2 (λ) along the path. Here G n,m (λ) = |E m (λ) − E n (λ)| is the energy difference between the m th and n th eigenstates and λ parametrizes the path, i.e. it represents the values of [Jκ] for the links on the path. The degeneracy decreases with the minimum distance d. The ratio G 2,1 (λ)/G 3,2 (λ) for small values of κ (i.e. from the set (i)) exhibits more involved behavior as the gap between the ground states and higher excited levels is much smaller in this case.
It is interesting to point out that fermionic parity of the system may change while the particles are exchanged. The average values of the parity of the system are shown in Figure 5(b) where 2 and 1 are assigned to the even and the odd parity sectors respectively. The parity of the system sets the small gap G 2,1 between the nearly degenerate states as G 2,1 = E 2 − E 1 for odd systems and G 2,1 = E 2 + E 1 for even systems, where E i denotes the spectrum of the system increasing with the index i ≥ 1. However, the parity does not affect G m,2 for any value of m. This is one of the causes of the oscillations seen in Figure 5. It is not the only cause though. Calculations show that the average energy of the nearly-zero modes (E 1 and E 2 ) oscillates in the same way as well.
For a system which is initially in the ground state space, the exchange of B and C transforms the system within the 2-dimensional ground state space provided that the process is adiabatic (Appendix A contains a detailed discussion on adiabaticity and the path which is followed). When the change of the parameters of the Hamiltonian follow a smooth closed curve in the parameter space, then the evolution of the system is governed by the following Berry phase (matrix) B(C) [30]: , 2}, and λ 1 and λ M +1 are coinciding points denoting the beginning and end point of the closed trajectory whose curve length length(λ M +1 , λ 1 ) is divided into M equal pieces ∆λ =length(λ M +1 , λ 1 )/M. Note that because the path traveled by vortices encloses zero area, all local and geometrical effects are eliminated and only the topological interaction is realized.
The Berry matrix B(C) is written in the same basis as that of A(λ 1 ) (see Appendix A for more details about λ 1 ). However, it is independent of the choice of the bases in which other A(λ)s are written, as long as the basis states are smooth functions of the parameter λ. For practical purposes, we used the nearly degenerate energy eigenstates (|Φ 1 (λ) and |Φ 2 (λ) ) of the system. We point out for clarity that |Φ 1 (λ) and |Φ 2 (λ are always distinguishable thanks to the small gap which separates them and which is never less than 10 −6 for all values of κ and d we used.

Numerical Methods for Calculating the non-Abelian Berry Phase
In this Section, we present the arguments that we use for the numerical calculation of the Berry matrix. The results and discussions are presented in the next Section.
In order to calculate the Berry matrix B(C) (20), we first need to evaluate There are many different ways to approximate the derivative of a function [31], but we are going to use the central-difference formula which says as long as the third derivative of f is continuous. It approximates f ′ (x) on the order O((∆x) 2 ); however it is also possible to get higher order approximations [31]. In our case, we performed the calculations of f ′ (x) on the order O((∆x) 6 ). By applying the formula (22) to A ab (λ k ), we get This reduces the Berry matrix calculation to finding the overlaps between the ground states of adjacent points on the trajectory. Let M be the number of data points used to calculate the Berry matrix (20). For all data points k = 1, ..., M and some other point k = 0, let |φ k be the ground states of β(k)-fermions defined as as in equation (14). Since the T 's are unitary matrices, we can write for k = 0 and replace c † ↔ c ↔ of Equation (23) with the left hand side of Equation (24) as where By applying the Bloch-Messiah theorem to Equation (25), we can only get the absolute value of the overlap between the ground states | φ 0 |φ k | = | det U(k, 0)| (see Appendix B). However, it is possible to get the complete overlap by using the Thouless theorem [28,32]. The Thouless theorem states that when φ k |φ 0 = 0, the ground state |φ k can be written as Having the ground state |φ k represented as |ψ (k,0) for k = 1, ..., M where we have fixed the overall phase to 1 (Equation (27)), we recall that the excited states can be represented by using the column exchange technique on T (k, 0) discussed at the end of Section 2. Then, the overlap between the ground states |ψ (l,0) and |ψ (k,0) for k = l = 0, reads In a recent paper [33], the overlap φ 0 |e Z † (l,0) e Z(k,0) |φ 0 has been calculated as where Pf denotes the Pfaffian, Z(l, 0; k, 0) = Z(k, 0) −I I −Z * (l, 0) and N is the number of fermions (and therefore also the size of Z). Since the numerical algorithms for calculating the Pfaffian of large matrices are generally slow, it is more convenient to proceed using the following expression (see Appendix C) Although the sign in Equation (28) is ambiguous due to the square root, for a series of smoothly varying matrices the correct sign can be traced from their previous values. We point out here that, for the purpose of calculating the Berry matrix (20), the degenerate ground states at each point of the trajectory in parameter space is defined in terms of the reference state |φ 0 . Thus the reference state must be chosen in such a way that it is not orthogonal to any of the ground states along the whole trajectory, as in Equation (27).

Numerical Results and Discussions
Before presenting the numerical results, let us briefly discuss the Ising anyons consisting of three different particles: 1 (i.e. vacuum), ǫ and σ. Bringing two particles together is called fusion, and Ising anyons satisfy the following fusion rules: where × denotes the fusion operator. For example, the first rule says that two σparticle may either annihilate or fuse into an ǫ-particle. Although the fusion of two Abelian anyons gives one outcome, non-Abelian anyons have multiple fusion possibilities or channels which are one way of accounting for the degeneracy of the ground state. Thus σ-particles are non-Abelian anyons while the other two particles are Abelian.
To detail the nontrivial implications of these rules, consider a system with four σ-particles (a, b, c and d) whose total topological charge corresponds to the vacuum. In this system, the fusion of any two σ-particles (say a and b) determines the fusion channel of the other two σ-particles (c and d) and because there are two different fusion results that a and b can fuse into, there is a two dimensional fusion space associated with the system. The basis of this space can be chosen as the resulting fusion states of a and b as {|Ω ab 1 , |Ω ab ǫ }. On the other hand, another basis can be chosen based on the fusion results of b and c {|Ω bc 1 , |Ω bc ǫ }. The matrix that relates the two fusion bases is called the F -matrix. It is analogous to the 6j symbols encountered in the couplings of three spin-1/2 particles. For arbitrary numbers of anyons, different sequences of fusion-basis transformations, starting from a particular fusion-basis ending in another one, must be equal. This imposes a consistency condition on F -matrices known as a pentagon equation [11].
On the other hand, particle exchange operators τ are represented by R-matrices on this fusion space. Nontrivial relations arise when we consider particle exchange operators together with F -matrices. These relations can be expressed by the hexagon equations [11]. Solving the pentagon and hexagon equations specify the consistent anyon models. For the Ising model, there are consistent solutions as follows R ǫǫ 1 = −1, R σǫ σ = R ǫσ σ = iα, R σσ 1 = θe iαπ/4 , R σσ ǫ = θe −iαπ/4 with the possible combinations of topological spin θ (associated with the counterclockwise rotation of a particle by angle 2π around itself) and α θ = e iπν/8 , α = (−1) (ν+1)/2 where ν is the spectral Chern number taking odd integer values [11]. For translationally invariant configurations of the Kitaev Honeycomb model (e.g. the reference X, Yconfiguration i.e. X q = 1 and Y q = 1 for all q), the spectral Chern number is equal to ±1 (depending on the direction of the magnetic field) [11].
The solution of the pentagon equation also determines the transformation relation between the fusion bases {|Ω ab 1 , |Ω ab ǫ } and {|Ω bc 1 , |Ω bc ǫ } of σ-particles (a, b, c and d) as up to an arbitrary relative phase e iϕ . In that regard, the representation R of the braiding operator which exchanges σ-particles b and c is diagonal for the fusion basis {|Ω bc 1 , |Ω bc ǫ } where {bc} stands for the basis {|Ω bc 1 , |Ω bc ǫ } and by using the transformation relations above, R {ab} reads On the other hand, the Berry matrix B(C) given by Equation (20) is written in a different basis. This is the same basis as A(λ 1 ), where λ 1 is the starting point of the trajectory C in the parameter space. These basis states consist of two nearly degenerate energy eigenstates {|Φ 1 (λ 1 ) , |Φ 2 (λ 1 ) } (see Sec. 3 also). We mention for completeness that the coupling coefficient of the starting configuration is slightly different than that of the original configuration shown in Figure 4. We discuss the details of the trajectory used in the calculations in Appendix A.
Because the Berry matrix is not invariant under the basis transformation and the relation between the fusion channels and the basis {|Φ 1 (λ 1 ) , |Φ 2 (λ 1 ) } is unknown, the comparison of the eigenvalues of the Berry matrix with {R σσ 1 , R σσ ǫ } is more meaningful. However, before we do that we would like to point out an interesting similarity between the calculated Berry matrices and an R-matrix: by making an analogy between four σ's (a, b, c and d) and four vortices (A,B,C and D) in Figure 4, for the values of κ from as the minimum distances d increases. We point out that the arbitrary relative phase e iϕ = −i is chosen for R {AB} ; and a similar arbitrary relative phase is fixed between the basis states |Φ 1 (λ 1 ) and |Φ 2 (λ 1 ) . This arbitrary phase reflects a gauge freedom and is chosen to obtain the best agreement between B(C) to R {AB} . The numerical results are summarized in Figure 6 where the Frobenius norm was used to measure the distance between the two matrices. The Frobenius norm of an n × n matrix A is defined as |A| F = i,j |A ik | 2 1/2 = tr(AA † ). Here, we calculate the norm |R {AB} − B(C)| F to measure the distance between B(C) and R {AB} . The high level precision of the results is reflected in the unitarity of the calculated Berry matrix: B(C)B † (C) − I F < 10 −5 for all d and κ. It is not difficult to show that the maximum Frobenius distance between two unitary 2×2 matrices is equal to 2, so to evaluate the error of the calculations the Frobenius distance has to be divided by 2.
As an example, for the vortex configuration illustrated in Figure 4  whose Frobenius distance to R {AB} given in Equation (5) is equal to 4.8 × 10 −5 and the normalized error associated with B(C) compared to the exact result from the effective field theory R {AB} is 0.0024%. This similarity between the basis {|Φ 1 (λ 1 ) , |Φ 2 (λ 1 ) } and {|Ω AB I , |Ω AB ǫ } can be understood from the perspective of Majorana fermions. First of all, note that the eigenstates |Φ 1 and |Φ 2 are also the eigenstates of the operators iγ 1a γ 1b and iγ 2a γ 2b where , are Majorana fermions, which are their own anti-particles. Kitaev showed that for odd values of the Chern number vortices carry Majorana fermions [11]. The localization of the Majorana fermions around vortices are also demonstrated numerically in [34] by expressing them in terms of c-fermions. In that respect, the equality of the bases can be understood as the Majorana fermion pairs (γ 1a , γ 1b ) and (γ 2a , γ 2b ) being localized on the vortex pairs (A, B) and (C, D). Note that because the localization of Majorana fermions is only related to the eigenstates of the Hamiltonian, it does not depend on the history of the vortex motion. That is, the configuration properties (such as relative distance of vortices to each other, boundary conditions, the values of L x and L y etc.) determine the relation between the energy basis and the fusion basis. For example, for the same size configurations with L x = 1 and L y = −1, the Berry phase B(C) is diagonal so that (Berry phase calculations for the other L x , L y and N x , N y configurations will be presented in a different context in future.) On the other hand, the eigenvalue comparison is presented in Figure 7 where The exponential convergence of the Berry phase to the expected values as the minimum distance d between the vortices increases, which the accuracy of our calculations allows us to observe, reveals subtle connections. The oscillations and the agreement of the calculated Berry phase to the expected values are in correlation with the exact ground state degeneracy and the gapfulness of the system ( Figure 5) with few exceptions. However, one should keep in mind that although the ground state degeneracy and the gapfulness of the system play a role in the physical evolution of the system, the Berry phase only depends on the eigenstates. Therefore, all the Berry phase values should be understood in terms of the properties of the eigenstates.
At this point, it is essential to note that the structure of the eigenstates determines the spread of the localization of the Majorana fermions around the vortices [34]. Because the braiding properties of the non-Abelian anyons are determined when the particles are away from each other to avoid tunneling, the radius of the Majorana fermions around the vortices affects the calculated Berry matrices. Thus, the relation between the calculated Berry phase to the expected R values is also an indication of how tightly the Majorana fermions are localized around the vortices.
In that respect, the correlations between Figure 5 and Figure 7 suggests that the bigger gaps and the nearer degenerate states make the localization of the Majorana fermions tighter around the vortices. However, it seems that this is not always the case as is seen for the systems with d = 6 and the small values of κ (i.e. from the set (i)).

Outlook and Summary
To summarize, we have analyzed some four-vortex configurations of the Kitaev honeycomb model (see Figure 4) under various magnetic fields. All lattice configurations are even-by-even in terms of the number of plaquettes and they are defined on a torus whose size is determined by the minimum distance d between the vortices. As was shown in Section 3, the ground states of such configurations are two-fold nearly degenerate. Under the adiabatic interchange of the vortices, the evolution of the system is restricted to the ground states and it is governed by the Berry matrix whose numerical evaluation was presented in detail in Section 4. The numerical results show that as d increases the small gap between nearly degenerate state decreases exponentially with some oscillations (see Figure 5) also the calculated Berry matrices get exponentially closer to the expected ones with similar oscillations (see Figure 7).
The main result of the paper is the explicit demonstration of the non-Abelian statistics by numerically calculating the Berry matrix which governs the evolution of the system under the adiabatic interchange of vortices of the Kitaev honeycomb model. We showed that the resulting Berry matrices exponentially converge to the expected braiding properties of Ising anyons derived from the effective field theory. The presented method for the Berry phase calculations represents a direct approach to the non-Abelian statistics, and it can be used to study braiding properties of anyons for any configurations. The high accuracy of the method also allows us to test the stability of the Berry phase, so that we may make conclusions about both its implementation under (more) realistic conditions and the pitfalls we should expect in such situations. Moreover, the accuracy of the calculations presented here can provide important predictive power in the analysis of other, less understood, topological phases. They also reveal the dependence on details of the implementation of the braiding operations and thus allow testing of non-adiabatic and other effects which will likely constitute a dominant source of errors in topological quantum information processing.
The accuracy of the calculations shows the exact dependence of the simulated Berry phase on the details of the model, like the splitting of the ground state levels intrinsic to any finite system. Possible applications thus extend to modeling and implementation of quantum information protocols whose reliability can be tested under various effects; for example, disorder. Naturally, the first step in this potentially fruitful story is the demonstration of highly accurate and sufficiently large scale direct simulation of the Berry phase as we have presented in this manuscript. Moreover, the method is general enough to be useful in the study of non-Abelian statistics in other models including the Yao-Kivelson model [21], square-octagon model [22] or any system with a quadratic fermionic Hamiltonian.  Figure A1(a). To move the vortex smoothly, we need to start changing [Jκ] y q+ nx before [Jκ] y q reaches −a, as shown in Figure A1(b). For this reason a vortex is never completely localized on a particular plaquette. We will refer to the part where only one [Jκ] changes as the linear part and the part where two [Jκ] changes as the circular part. In our numerical calculations, every link on the trajectory is sampled with 4000 data points (3991 linear plus 9 circular) which are separated by an equal distance ∆s in the parameter space. Let l be the length of the linear part from 0 to the beginning of the circular part (see Figure A1(b)) and r the radius of the circular part; then r + l = |a| and ∆s = ∆l = r∆θ where ∆l = l/3991 and ∆θ = π 2 /9, so that l = 0.997|a| and r = 0.003|a|.
The starting point of the trajectory (i.e. λ 1 ) is chosen to be the beginning of the linear part of the [Jκ] of the first link (i.e. [Jκ] y (2,Ny ) ) of the path; because of that the coupling coefficients of the starting configuration are slightly different than that of the original configuration shown in Figure 4.
On the other hand, the first requirement for a process to be adiabatic is the gap between the eigenspace and the rest of the Hamiltonian's spectrum and this has been analyzed in Figure 5. The standard condition for maintaining the adiabaticity for the n th eigenstate is as follows: where dot denotes the time derivative. However, this condition is only valid if Φ n (λ)|Φ m (λ) and E m (λ) − E n (λ) are constant for all m through the trajectory [29]. When Φ n (λ(t))|Φ m (λ(t)) and E m (λ) − E n (λ) are not constant the validity condition where v = dλ dt is the speed of the process and we used the following equality Φ n (λ(t))|Φ m (λ(t)) = − Φ n (λ(t))|Φ m (λ(t)) which can be verified by differentiating the orthogonality condition Φ n (λ)|Φ m (λ) = 0. First we define for a k-fold degenerate eigenspace with basis |Φ n i (λ) for i = 1, ..., k. Then by using the following inequality we can rewrite the general validity condition for adiabaticity as v max d dλ Φ n i (λ) ⊥ min |E n (λ) − E m (λ)| ≪ 1. Figure A2 shows the validity of the adiabatic approximation for various κ values against d. The smaller the values of κ, the smaller the gap to the first excited state, and thus the smaller the speed v is required to maintain adiabaticity. and express the overlap between |ψ (k,0) and |ψ (l,0) as ψ (l,0) |ψ (k,0) = | det U(l, 0)| | det U(k, 0)| (−1) N (N +1)/2 det U(l, k) det U † (k, 0) det U(l, 0) = (−1) N (N +1)/2 exp {iθ 0 (l, k)} | det U (j, i) |, where θ 0 (l, k) = arg det U(k, 0) det U † (l, 0) det U(l, k) .