The Kitaev honeycomb model on surfaces of genus $g \geq 2$

We present a construction of the Kitaev honeycomb lattice model on an arbitrary higher genus surface. We first generalize the exact solution of the model based on the Jordan-Wigner fermionization to a surface with genus $g = 2$, and then use this as a basic module to extend the solution to lattices of arbitrary genus. We demonstrate our method by calculating the ground states of the model in both the Abelian doubled $\mathbb{Z}_2$ phase and the non-Abelian Ising topological phase on lattices with the genus up to $g = 6$. We verify the expected ground state degeneracy of the system in both topological phases and further illuminate the role of fermionic parity in the Abelian phase.


Introduction
The Kitaev honeycomb model is an example of an exactly solvable two-dimensional model that exhibits both Abelian and non-Abelian topological phases [1]. The Abelian phase, which is also known as the toric code [2], provides a realization of a topological quantum field theory known as doubled- 2 theory. The non-Abelian phase is effectively described by the Ising topological quantum field theory [3]. The main attribute of topological field theories is a dependence of the dimension of the relevant Hilbert space on the topology of the underlying manifold on which these theories are realized. For example the doubled- 2 theory is represented in the twodimensional toric code by a non-degenerate ground state on a genus 0 surface like an infinite plane or a sphere, and a four-fold degenerate ground state on a genus 1 surface like a torus. Similarly, the Ising topological field theory is linked to a three-fold degenerate ground state of the honeycomb lattice model in its non-Abelian phase on a torus. However, the square lattice of the toric code and the honeycomb lattice of the Kitaev model permit realizations of only these two surface topologies as the Euler characteristics of both lattices are zero.
We extend the solution of the Kitaev honeycomb model to closed surfaces of genus greater than one. We will rely on the exact solution of the model based on Jordan-Wigner fermionization [4][5][6]. This solution allows us to factorize the model into a fermionic superconductor on a topological doubled- 2 square lattice background or vacuum state. In order to generalize this to higher genus surfaces, we introduce a lattice which can be realized on such surfaces and accordingly adjust any definitions and relations to this context. We will then demonstrate the generalized solution on a number of different surfaces of genus greater than 1 by calculating the ground state degeneracy of the model in both the Abelian and non-Abelian phases. In this context, we also investigate additional features of these topological states that are intrinsic to their lattice realizations.
A natural framework for our investigation of two-dimensional lattice models whose topological phases effectively realize certain topological quantum field theories is the axiomatic definition of these theories. n-dimensional topological quantum field theory is defined as a functor from a category of n-cobordisms to a category of vector spaces [7,8]  -( ) F n : Cob Vect. 1 is subject to certain axioms which for example ensure that vector spaces originating from topologically equivalent manifolds are isomorphic and that the disjoint union of (n−1)-manifold carries over to a tensor product between vector spaces. The functor satisfying these axioms is called modular and the underlying categories are called monoidal. We point out that realizing a topological phase of a physical system on a closed oriented surface of some genus represents a realization of an important part of this functor. Specifically it assigns to the surface (2-manifold) a vector space spanned by the ground states of the relevant physical system. We first give a concise overview of the model and its effective spin/hardcore-boson representation on a square lattice in section 2. A realization of lattices on higher genus surfaces is introduced in section 3, which is then followed by the implementation of the model on these lattices and its solution using the Jordan-Wigner fermionization in its effective spin/hardcore-boson representation in section 4. The last two sections describe the calculation of the ground state in section 5 and evaluation of the ground state degeneracy of the model on surfaces with the genus 2-6 in section 6.

The model
The Kitaev model [1] is a honeycomb lattice with a spin 1 2 particle attached to each vertex. Each spin interacts only with its nearest neighbors via an interaction term that depends on the orientation of the link (x, y or z) connecting them. Explicitly, if i and j label neighboring vertices connected by a link of orientation α, these spins interact via a term of the form s s We can also add a time-reversal and parity-breaking potential to this Hamiltonian which comes from third order perturbation theory of a weak magnetic field. The effective potential is k = å V V p p where κ is a coupling constant and the sum is over the plaquettes of the system with each hexagonal plaquette making the following contribution to the potential: where the sites of the plaquette p have been numbered as in figure 1. Hence, the full Hamiltonian of the model is We can define a vortex operator W p for each plaquette p of the lattice. If we number the sites of the plaquette p as in figure 1 The vortex operators W p commute mutually and with the full Hamiltonian, including the time-reversal and parity-breaking potential terms. Consequently the Hilbert space can be written as w p is the common eigenspace of the W p operators corresponding to the particular configuration of eigenvalues {w p } where w p =±1. We say that a vortex occupies the plaquette p if w p =−1 [1,9].
Kitaev [1] solved the system by a reduction to free fermions in a static  2 gauge field. He showed that the model exhibits four distinct topological phases including three Abelian toric code phases A x , A y , A z , satisfying Figure 1. The Kitaev honeycomb lattice model and its phase diagram. As shown on the left, any given link has one of three possible orientations, x, y or z and we number the vertices of a plaquette 1-6. The phase diagram can be thought of as the convex hull of the three points (J x , J y , J z )=(1, 0, 0), (0, 1, 0), (0, 0, 1).

:
, : , : , 6 x y z x y x z y z x y z and an additional phase B which occurs when all three inequalities above are not satisfied simultaneously. In the absence of the magnetic field the B phase is gapless, but in the presence of a magnetic field it acquires a gap and becomes the non-Abelian Ising phase. Its quasi-particle excitations, which in our representation are formed by Majorana fermions attached to vortices, show non-Abelian fractional statistics and are known as Ising anyons. As described in [6], the model can be mapped onto a square lattice whose vertices carry effective spins and hardcore bosons. In this representation, the Hamiltonian of the model (2) acquires the following form x q q q q n x q n q nx y q q z q q q n y q n q ny where the t a q are the Pauli operators for the effective spin at a site q and † b q and b q are creation and annihilation operators for hardcore bosons. The sums in the Hamiltonian are over all the sites of the lattice. The contribution to the potential V of a plaquette P has the following form in this representation The vortex operator W P for each plaquette P of the lattice is now defined as t t t t = -- q q q is the boson number operator. We say P is occupied by a vortex if the eigenvalue of the corresponding vortex operator is −1 and is empty otherwise.

Lattices on higher genus surfaces
We will now discuss the construction of lattices on closed surfaces with different topologies. To define the model on a closed surface of a higher genus, g>1, we necessarily have to consider a lattice with an Euler characteristic χ that is negative due to the relation c =g 2 2 . A perfect square lattice with the number of vertices V=N, the number of plaquettes or faces F=N and the number of edges = E N 2 permits at most a closed surface of genus g=1 as its Euler characteristic χ=V−E+F is zero. To construct a lattice with negative Euler characteristic from a square lattice requires alterations of some of its vertices or plaquettes. For example, they may include increasing or decreasing the number of edges connected to some vertices or changing the number of edges associated with some plaquettes. We refer to such alterations as defects and emphasize that these are local lattice defects as opposed to non-local defects such as lines of dislocations [10,11,12]. We can think of these defects as particles, called genons [13,14].
We first construct a lattice with g=2 before considering lattices of higher genus. We start with an octagonal piece of square lattice and identify or glue its opposing boundaries together in a way similar to creating a torus by identifying opposite boundaries of a rectangle. The construction is illustrated in figure 2. If we tessellate an octagon with a square lattice and identify the sites residing on the boundary as indicated in figure 3, the resultant lattice will be of genus g=2. We now have a defect plaquette with 12 edges centered around the corners of the original octagon, which are all identified once the boundary edges are glued together. Clearly we could tessellate an octagon with a variety of square lattices of different sizes. The particular lattice we use is characterized by three numbers {N a , N b , N c } which specify the number of vertices along the vertical, diagonal and horizontal edges respectively as shown in figure 3. The total number of vertices on such a lattice with dimensions c a c tot . We can calculate the Euler characteristic by noticing that there are exactly 2N tot edges and N tot −2 plaquettes including the defect plaquette. Hence we have c = - as desired. We note for completeness that there are other ways of gluing the edges of an octagon together in order to produce a g=2 surface but these may lead to the emergence of undesired line defects. Our approach described above avoids this issue.
Alternatively, we could consider a similar construction using the dual lattice. While the original lattice has each vertex four-valent and all plaquettes are square except the defect plaquette, the dual lattice has all plaquettes square and all vertices four-valent except one defect vertex which is twelve-valent. However, this would require changes of the Hamiltonian of the model. We therefore prefer to work with the original lattice which preserves the form of the Hamiltonian. We will, nevertheless, need to define the vortex operator for the defect plaquette and its magnetic contribution.
We now consider the construction of lattices on surfaces with genus g>2. One approach to generalize the construction developed for the g=2 surfaces above would be to start with a polygon with a greater number of sides (e.g. dodecagon for a g=3 surface) and then glue the opposite sides accordingly. Here we prefer a different and more modular approach which lends itself more naturally to a numerical implementation.
Consider the octagonal piece of lattice as described above. Once all but two of the edges have been glued together, we are left with a lattice with the topology of a torus with two punctures in it. We now use this as a building block for constructing lattices with higher genus. Consider g−1 copies of a torus with two punctures. We can always glue the punctures together in such a way that results in a closed connected surface of genus g. With regards to the lattice, we start with g−1 copies of the octagonal piece of square lattice described above and stitch them together to form a chain of octagons as depicted in figure 4. We now form a lattice on a surface of the desired topology by identifying the remaining opposing edges of each octagon as well as by gluing together the remaining edges of the first and the last octagon of the chain. The resultant lattice will have g−1 defect plaquettes, identical to the one described above, located where two octagons are joined together.
We now verify the Euler characteristic for our higher genus lattices. If each octagonal piece has dimensions N a , N b and N c then the total number of vertices on this lattice is (g−1)N tot . We can still uniquely associate every vertex to two edges so the lattice has 2(g−1)N tot edges. To write down the number of plaquettes as a function of the lattice  The lattices we will be considering on genus 2 surfaces will tessellate an octagon as depicted on the left. They are characterized by three numbers N N , a b and N c . N a is the number of links crossing the vertical (green) edge of the octagon. N b is the number of sites living on a diagonal (blue or red) edge. N c is the number of links crossing the horizontal (purple) edge of the octagon. When the edges have been identified appropriately, the links colored red in the center image form a closed chain depicted in the image on the right. The corners of the octagon all meet at a common point represented by the black dot at the center of the right image. The corresponding plaquette, centered around this point, will have 12 edges and we will refer to it as the defect plaquette. dimensions {N a , N b , N c }, we can associate every vertex to the upper right hand plaquette it forms a corner of. Every square plaquette will be assigned a unique vertex while the defect plaquettes will be assigned three. So the number of plaquettes on the lattice is equal to the number of vertices minus 2 for every defect:

The model on surfaces of genus g2
We now consider the model on the lattices constructed in the last section. We first write down and discuss the Hamiltonian for the system and its symmetries in the effective spin/hardcore-boson representation of the model. We then fermionize the bosons to obtain a Hamiltonian quadratic in fermionic operators. Since the lattices we will be considering do not have any translational symmetries, we will not be able to write down the ground state in closed form as was done in [6] for the model on a torus. However, the formalism allows one to efficiently diagonalise the Hamiltonian numerically within any particular common eigen-subspace of the models symmetries.
The Hamiltonian in the effective spin/hardcore-boson representation on the lattice described above is of the same form as that on a lattice without defects (2). In both cases, every vertex is four-valent with two horizontal (x-links) and two vertical (y-links) edges attached. If we denote a site of the lattice by q, then by q+n x we denote the neighbor to the right of q that is connected to it by an x-link. Similarly, we use the notation q+n y to denote the neighbor above q that is connected to it by a y-link. The bare Hamiltonian can then be written as follows: Regarding the potential k = å V V P P , the contribution from the square plaquettes are still given by the expression (8). On the other hand, the contribution of the defect plaquettes to the potential are more complicated and are given by This expression follows from translating the three-body spin terms, linked at the third order of perturbation theory to the weak magnetic field, into the effective spin/hardcore-boson representation. In the original honeycomb picture, the defect corresponds to a plaquette with eighteen edges and the contribution to the potential by a defect plaquette is the following sum of three-body spin terms: where the sites of the plaquette are numbered as depicted in figure 5.
We can still define a vortex operator which commutes with the full Hamiltonian k = + å H H V P P 0 for every plaquette. For square plaquettes the vortex operator is defined as in equation (9). For the defect plaquettes however we define the vortex operator as follows t In addition to the vortex operators, we also define an operator, which commutes with the Hamiltonian, for every generator in a basis for the 1st  2 -homology group H 1 of the lattice. We will call these operators loop operators and to define them we will need to choose a basis for H 1 and a particular representative from each homology class in that basis. For a lattice of genus g2 the rank of H 1 is g 2 so we will need to choose g 2 homologically distinct cycles. We will choose the cycles depicted in figure 6 and their associated homology classes as the representatives and basis respectively. As depicted, for a lattice with g−1 copies of an octagonal piece of lattice we will choose three cycles on the first copy, two cycles on every other copy (reflecting the fact that every additional copy increases the genus by 1 and the rank of H 1 by 2) and one horizontal cycle that spans each octagon. The loop operators we will define for these cycles will act on the sites of the lattice that are connected to the links that constitute the cycles. How a loop operator acts on a particular site is determined by the way the associated cycle passes through it. There are six ways a cycle can pass through a site as depicted in figure 7 and we will associate a single site operator with each of them as follows: . 15  vertical part and two corners and so the loop operator for this cycle can be written as follows:

t t t t t t t t t t t
Corner 1

Horizontal
Corner 4

Vertical
One can in principle define a different set of loop operators that commute with the Hamiltonian but these will in general be equivalent to a product of the loop operators already defined times a product of vortex operators [9]. The vortex and loop operators form a set of commuting observables, allowing us to decompose the Hilbert space as follows:   subspaces where it can be expressed as a combination of terms that are quadratic in fermionic operators. The restricted Hamiltonian can then be diagonalized by an appropriate Bogoliubov transformation. We now change to a basis of the Hilbert space which reflects the decomposition (17). It seems natural to consider the common eigenvectors of the vortex and loop operators along with the eigenstates of the boson number operator ( = † N b b q q q ) for each site q of the square lattice. However, the basis so defined would be overcomplete. If there are -( ) g N 1 tot sites in the lattice, the model clearly has distinct combinations of eigenvalues a common eigenvector of this set of observables might have. However, the vortex and number operators are not completely independent operators as they satisfy two conditions.
The first condition is the fact that the product of all vortex operators is equivalent to the identity operator. That is, where the product is over all the plaquettes of the lattice. Since a product of vortex operators can be thought of as counting the parity of vortices occupying the associated plaquettes, this essentially means there can only be an even number of vortices in total in the model. So the number of independent vortex operators is and hence the number of configurations of vortices in the model is - The second condition is a relation between the parity of bosons in the system and a certain product of vortex operators. For a lattice where the numbers N a and N b are both even, we can consider a set of plaquettes forming a checker board pattern as depicted in the top left image of figure 8 by the colored squares. It is easy to check that since the Pauli operators square to the identity, the product of the vortex operators associated with the colored (or uncolored) plaquettes is equivalent to the boson parity operator where q runs over all the sites of the lattice. In other words, the parity of the number of bosons must be the same as the parity of the number of vortices on colored plaquettes (or equivalently uncolored plaquettes). Since the parity of bosons is fixed to be 1 or −1 depending on the configuration of vortices, the number of independent boson number operators N q is -- tot and hence the number of configurations of bosons in the model is - tot . For lattices where N a or N b are odd numbers, there is a similar dependence of the boson parity on the configuration of vortices in the system. For such lattices, we cannot color the plaquettes with a perfect checker board pattern but we can consider sets of plaquettes as depicted in figure 8, such that the checker board pattern is misaligned along a 1-cycle of links that separate plaquettes of the same color. The exact pattern we choose for coloring in plaquettes and the associated cycle along which the checker board pattern is misaligned depends on the parity of the numbers N a , N b and g for the lattice and is described in figure 8. If we compose the corresponding vortex operators, the Pauli operators for sites away from this cycle will cancel out as they did before but along the cycle, the resultant operator will act with a string of Pauli operators and may not act with the parity operator -( ) N 1 2 q for some sites. However, we can cancel the action of these Pauli operators, and replace any missing single site parity operators we need to obtain the full boson parity operator, by composing this product of vortex operators with a product of loop operators that act on the sites connected to the links of the cycle. The desired product of loop operators that act on the sites connected to the links of the cycle are shown in figure 8. In general, the boson parity operator can be written as The next step of the solution is to use a 'Jordan-Wigner' type transformation to fermionize the bosons of the model. This should result in a Hamiltonian which is quadratic in fermionic operators which we will then be able to solve using the Bogoliubov-de Gennes (BdG) technique. To fermionize the bosons, we will define a Jordan-Wigner type string operator S q for each site q of the lattice. The composition of these string operators with the boson creation and annihilation operators will be fermionic creation and annihilation operators. Expressing the Hamiltonian and other observables in terms of these new operators will effectively transform the hardcore bosons of the model into fermions.
To define a string operator for a site q of the lattice we consider the following: if we had a particle located at the reference site (as in figure 9(a)) we can always move that particle to any site q by first moving it to the right an appropriate number of sites and then up an appropriate number of sites. Even the sites below the level of the reference site can be reached in this way by making use of the boundary conditions as shown in figure 9(b). We can associate a single site operator for every site traversed in the path just described connecting the reference site to the site q. The string operator for q, denoted S q , will be defined as the composition of these operators. To every site i crossed by the horizontal part of the path we associate the operator t --( ) N 1 i i x , to the corner of the path we associate the operator t i y , to every site i crossed by the vertical part of the path we associate the operator t i x and to the last site of the path we associate the operator t i y . Since each of these operators act on different sites, they all commute with each other and so we are free to define S q as the composition of these operators without worrying about the order of composition. If we let q x denote the number of sites that need to be traversed in the horizontal part of the path with the site at the corner and q y the number of sites that need to be traversed in the vertical part of the path with the site at the end, then we can number the sites of the path from 1 to q x +q y , beginning at the reference site and ending at q and we can write the string operator for q as follows: If we consider two string operators S q and ¢ S q such that ¹ ¢ q q there will be a single site, shared by the paths defining the string operators, where the action of S q anti-commutes with the action of ¢ S q . It follows that composing the string operator S q with the bosonic creation and annihilation operators for the site q defines fermionic creation and annihilation operators for q which we denote by † c q and c q .
, , 0, , 0. 23 Expressing the basic Hamiltonian in terms of these fermionic creation and annihilation operators yields the following sum of quadratic fermionic terms, Noting that both the ¢ X q q , and ¢ Y q q , operators, being products of string operators, act on a closed loop of sites, we can associate a 1-cycle with each of the ¢ X q q , and ¢ Y q q , operators, namely the set of links joining the sites being acted on. These operators will always be equivalent to a product of loop operators, which is determined by the homology class of this cycle, and a product of vortex operators which is determined by a certain 2-chain related to the homology class of the cycle and the representatives of the homology classes we have chosen as a basis for H 1 . Recall that each loop operator is associated to a non-trivial cycle, the collection of which represent the generators of H 1 . So whatever the homology class may be for the cycle a associated with an X or Y operator, we can always create a unique cycle b which will be homologous to a by adding some combination of the cycles associated with the loop operators. A particular X or Y operator is proportional to the product of the loop operators corresponding to the cycles used in the combination forming b.
There will also be a 2-chain, which we denote by ς, which will have a+b as a boundary. A particular X or Y operator is also proportional to a product of the vortex operators associated with the plaquettes which constitute ς. We note that while such a 2-chain ς is not unique, the operator obtained by multiplying the vortex operators associated with the plaquettes of the 2-chain ς is unique. For example, if we cut out a cylinder with boundaries a and b from a torus, the product of the vortex operators inside of the cylinder is the same as in its complement. This follows from the fact that vortex operators square to the identity and the relation (18). In general, when expressed in terms of loop and vortex operators, the X and Y operators are of the same form. We will use a notation to reflect this by letting Z q denote X q if q is a x-link and Y q if q is a y-link. Explicitly we have 1 is the homology class of the cycle associated with the link q described above. So the Hamiltonian can be written as follows  tot tot matrices ξ and Δ are given by Here, d ¢ Regarding the potential, when expressed in terms of the fermionic creation and annihilation operators, each term appearing in the sum defining the contribution from a plaquette inherits a product of string operators similar to the X and Y operators. The potential also becomes quadratic in fermionic operators and can be written as To be able to calculate the ground state energy numerically for a particular vortex/homology sector, we need to understand how the matrix T represents the state fñ | and how it can tell us the parity of the of the number of occupied c-fermions modes it has. In general, T will be of the following form: tot tot matrices which, since T must be unitary, must satisfy Bloch and Messiah were able to show that a unitary matrix of the form (41) can be decomposed as follows [15,16] where the -´- tot tot matrices D and C are unitary and bothŪ andV are real matrices of the following block diagonal form: odd in half of the homology sectors and even in the other half. This leads to a splitting in the energy between fermionic ground states in half of the homology sectors from the other half resulting in the degree of degeneracy d=8. In figure 10(b), we see that the system with odd dimensions N a and N c has half of its homology sectors forming the ground state in the Abelian phase while the other half are excited states. As the system approaches the phase transition, the sectors forming the ground state begin to split with two of them becoming excited in the non-Abelian phase while four of the excited sectors drop in energy to join the remaining six non-excited sectors to form the ten-fold degenerate ground state in the non-Abelian phase. Due to finite size effects, there is a small splitting in the energy between the degenerate homology sectors that form the ground state. We expect this spitting to vanish in the thermodynamic limit. We measure this splitting by the difference in energy between the sector with the highest energy and the sector with the lowest energy. In figure 11 we plot the splitting between the the degenerate states as a function of N=N a =N b =N c for the two Abelian cases (even and odd sizes) and the non-Abelian case. As shown in the figure, we find the splitting between the sectors forming the ground state approaches zero exponentially as N grows. This calculation was done with κ=0.2 in each case and with J=0.1 for both of the Abelian cases and J=1 for the non-Abelian case. Figure 10. In (a) we show the difference between E min and the energy E of fermionic ground states and first excited states in each homology sector as a function of J x =J y =J for N a =N b =N c =4 and κ=0.2 on a genus 2 lattice. In (b) the same energy difference is plotted for N a =N b =N c =5. The number of degenerate ground states is included just above the lowest curves in both the Abelian and non-Abelian phases. Figure 11. The splitting in energy between the degenerate homology sectors, measured by the difference between the sector with the highest energy and the sector with lowest energy, vanishes exponentially as the system size N=N a =N c increases. The largest system size corresponds to over 5000 spins of the original honeycomb lattice; beyond that the numerical precision starts competing with the ground state splitting.
We used this method to calculate the degeneracy of the system on lattices with genus g=2, 3, 4, 5 and 6 in both the Abelian (for even and odd sizes) and non-Abelian phases. We have summarized the results in table 1. Kitaev showed using perturbation theory that the honeycomb model in the Abelian phase is equivalent to the toric code which can be shown to have a ground state degeneracy of 4 g . This agrees with our results for systems with even N a and N c . For systems where either N a or N c are odd we find the degeneracy is exactly half of 4 g . This can be attributed to the fact that the equivalent toric code in this case has a line defect in it like the one discussed in [18]. Our results for the non-Abelian phase agree with a formula discussed by Oshikawa et al [19] who showed that the bosonic Pfaffian state, which belongs to the same universality class, has a ground state degeneracy + -( ) 2 2 1 g g 1 given by the number of even spin structures on a surface of genus g [20].

Conclusion
In summary, we realized the Kitaev honeycomb model on surfaces with genus g2 by introducing extrinsic defects to the underlying lattice. This required a non-trivial generalization of the exact solution of the model to include extra loop symmetries associated with homologically non-trivial loops which are introduced by increasing the genus of the lattice. We also highlight the dependence of the fermion parity on both the vortex and loop symmetries of the model for various lattice dimensions. The generalized solution was then used to calculate the ground states in both the Abelian and non-Abelian phases of the model. The degree of degeneracy of these ground states in both topological phases are in accord with available theoretical predictions based on topological quantum field theory.
Our work provides a direct realization of two distinct topological quantum field theories, specifically the Abelian doubled- 2 and non-Abelian Ising theory, on closed surfaces of higher genus. As such it provides a solid basis for further investigation of the model on various manifolds, including also manifolds with boundaries which would extend previous studies of the Kitaev model [21]. Recent works on time-dependent simulation of creation and annihilation of vortex-like excitation on defects in the Kitaev model on torus [10] suggest the possibility of a dynamical process where creation and annihilation of extrinsic defects would result in dynamical change of the model genus and thus its topology. Interestingly this incarnation of topological field theory would be close to its axiomatic definition as a modular functor from a monoidal category of cobordisms to that of vector spaces [4,5].