Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model

We study the collective states of interacting non-Abelian anyons that emerge in Kitaev's honeycomb lattice model. Vortex-vortex interactions are shown to lead to the lifting of the topological degeneracy and the energy is discovered to exhibit oscillations that are consistent with Majorana fermions being localized at vortex cores. We show how to construct states corresponding to the fusion channel degrees of freedom and obtain the energy gaps characterizing the stability of the topological low energy spectrum. To study the collective behavior of many vortices, we introduce an effective lattice model of Majorana fermions. We find necessary conditions for it to approximate the spectrum of the honeycomb lattice model and show that bi-partite interactions are responsible for the degeneracy lifting also in many vortex systems.


Introduction
One of the main trends in contemporary condensed matter physics is the study of topological order. Apart from their fundamental interest, the motivation derives from the discovery of topological quantum computation [1,2]. The quasiparticle excitations of topologically ordered systems are anyons that can be employed to perform fault resilient quantum information processing [3,4]. It is of interest to find out what are the sufficient conditions for topological order to exist, which systems are useful for quantum computation and how anyons emerge. While the research concentrated initially on the experimentally challenging fractional quantum Hall states or abstract spin lattice constructions, the recently discovered topological insulators have given rise to novel constructions of anyon supporting systems [5]. The experimental accessability of the latter has led to the theory and detection of Majorana fermions becoming the subject of intense research.
Majorana fermions are the simplest possible non-Abelian anyons [6]. While they are not universal for quantum computation [7,8], their realization and the subsequent demonstration of topological information processing is considered a significant stepping stone towards a full scale implementation. Physical systems that are predicted to give rise to Majorana fermions include p-wave superconductors [9], Kitaev's honeycomb lattice model [10] and its generalizations [11,12,13,14,15], the Moore-Read state in fractional quantum Hall systems [16] and lately also various superconductor/insulator interfaces [17,18,19]. Experiments have been performed on fractional quantum Hall states [20] as well as on topological insulator interfaces [5], but the existence of Majorana fermions remains still ambiguous. For the detection of Majorana fermions, and for their future applications, it is crucial to understand how their anyonic properties are manifest in the microscopic physics. The essential property to all topological quantum computing schemes is the existence of non-local topological degrees of freedom. For Majorana fermions they appear as fusion degrees of freedom that describe the different locally indistinguishable ways the anyons can behave collectively. However, these states can become locally distinguishable due to tunneling processes that lead to the degeneracy being lifted [21]. This has been studied in Kitaev's honeycomb lattice model [22], in the Moore-Read state [23] as well as in p-wave superconductors [24,25,26]. It has also recently been discovered that in many anyon systems the tunneling processes can not only lift the degeneracy, but can even drive transitions between topological phases [27,28,29]. The physical conditions under which this occurs, however, depend on the microscopics of a specific model. Only when the tunneling is exponentially suppressed, the topological phase is stable and its anyonic excitations can be employed for topologically protected quantum computing.
In the first part of the paper we analyze the degeneracy lifting in the non-Abelian phase of Kitaev's honeycomb lattice model [10]. By simulating continuous vortex transport [30], we extend the vortex-vortex interaction analysis of [22]. We find the anticipated oscillations in the energy splitting [24], show how they depend on the microscopics of the model and characterize the stability of the topological low energy spectrum by obtaining the relevant energy gaps. In the second part we introduce an effective Majorana fermion lattice model, an extension of the tight-binding model considered in [31], to study the degeneracy lifting in many vortex systems. Simultaneous interactions between vortices lead in general to a hybridization of the states corresponding to fusion degrees of freedom [32]. We find necessary conditions, i.e. relevant effective flux sectors, for the Majorana model to approximate this subspectrum. The flux sectors are interpreted to arise from the vortices of the underlying model, which provides a direct connection between the microscopic and the effective models. Furthermore, we find that the effective tunneling amplitudes are well approximated by the energy splitting due to vortex-vortex interactions. This confirms that the picture of freely tunneling Majorana fermions applies also in many vortex systems.
This paper is organized as follows. In Section II we review the solution of the honeycomb lattice model using Majorana fermionization. We show how it gives rise to an equivalence between coupling and vortex configurations and how this can be employed to simulate continuous vortex transport. In Section III we analyze the oscillating degeneracy lifting of the fusion degrees of freedom when the separation between two vortices is varied. We obtain the physical parameter dependence of both the oscillations and the energy gaps characterizing the topological low energy spectrum. In Section IV we introduce a lattice model of Majorana fermions and employ it to study the collective behavior of up to eight vortices.

The honeycomb lattice model
In this section we briefly review Kitaev's honeycomb lattice model and its solution by mapping it to free Majorana fermions. For a more rigorous treatment we refer to the original work [10] and to the subsequent developments [22,33,34,35,36,37,38]. We show how to relate the manipulation of the vortices to the manipulation of the model's physical parameters. This enables the simulation of continuous vortex transport [30], which is later used to study the physics of the anyonic vortices.

Solution by mapping to free Majorana fermions
Kitaev's model [10], 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 real nearest neighbour couplings on links of type α specified by their orientation (see Figure 1(a)). The second term is an effective magnetic field of magnitude K, which explicitly breaks time-reversal invariance. The sum in this term runs over all next to nearest neighbour triplets as described in [22]. The model can be mapped to a free fermion problem when the spins are represented by Majorana fermions [10]. This procedure reduces the Hamiltonian to the quadratic form Here c i = c † i are Majorana fermions that satisfy {c i , c j } = 2δ ij . The first term describes their hopping between nearest neighbour sites while the second term describes next to nearest neighbour hopping between sites shown in Figure 1(c). The Hermitian operatorŝ u ij =û † ij live on the links of the honeycomb lattice and should be understood as a Z 2 gauge field. Since [û ij , H] = 0, they have no dynamics. The gauge field interpretation follows from the local constraint D i | Ψ = | Ψ that the eigenstates | Ψ of the original Hamiltonian (1) have to satisfy for all sites i [10]. The operator D i anti-commutes with the three link operators acting on vertex i and thus it can be interpreted as a local gauge transformation. The Hamiltonian (2) is gauge invariant.
At the heart of the exact solvability of the model are the local symmetries [ŵ p , H] = 0, where, as illustrated in Figure 1 are Hermitian plaquette operators associated with every plaquette p. In the fermionized picture these becomeŵ that in the gauge theory language can be understood as gauge invariant Wilson loop operators. Their eigenvalues w p = −1 are interpreted as having a π-flux vortex on plaquette p. The pattern of all plaquette operator eigenvalues w ≡ {w p } is referred to as a particular vortex sector, that is created by fixing the gauge, i.e. choosing the pattern of the link operator eigenvalues u ≡ {u ij }. The vortex sectors label the physical sectors of the model. We consider finite systems of 2L x L y spins with L x and L y plaquettes in directions n x and n y (see Figure 1(d)). Due to particle-hole symmetry in the problem, diagonalization of (2) in vortex sector w yields in general the double spectrum for some complex coefficients [α w i ] j , are complex fermionic modes. When the diagonalization is performed numerically, ±E w i and α w i follow from solving the eigenvalue problem ih w α w i = −E w i α w i . Here h w is the real kernel matrix (2) when restricted to vortex sector w.

Gauge/coupling configuration equivalence and vortex transport
Our aim is to study the physics of the vortices by considering the spectral evolution as a function of the vortex sectors w. While the sectors form a discrete set, we can effectively interpolate between them, i.e. simulate the transport the vortices, as follows. As can be seen from (2), the local value u ij of the gauge field appears always uniquely paired with a local coupling J ij . This implies that while the vortex sector depends on the gauge through the plaquette operators, (3), from the point of view of the Hamiltonian, (2), we can regard u as the sign configuration of the couplings. It follows that the Hamiltonians H w and H w ′ on two different vortex sectors can always be connected by flipping the signs of some of the couplings, i.e. one can always find some coupling configurations J and J ′ such that H w (J) = H w ′ (J ′ ) ‡. Therefore, by manipulating the local couplings we can interpolate the spectrum continuously between different vortex sectors and thereby study the spectral evolution when vortices are created, transported and annihilated. When this is performed adiabatically, the vortex sector w does not strictly speaking change, but in our simulation the spectrum will evolve as if it had. This same method has been employed to evaluate also the non-Abelian statistics of the vortices [30].
From now on we adopt a perspective where the system has been initially prepared in the ground state belonging to the vortex-free sector [39]. The vortex sectors w are viewed as coupling configurations J ≡ {J ij } of equal amplitude, but with locally varying signs. To transport a vortex between two plaquettes or to create/annihilate a pair of vortices on them, one needs then to change the sign of the physical coupling J ij on the shared link. This can be performed continuously by tuning J ij → 0 → −J ij §, which we simulate by changing the sign of the coupling in S steps such that at step s its value is given by J ij (1 − 2s S ). As S becomes large, this approximates well a continuous process. The resulting transport protocol is illustrated in Figures 2(a)-(d). Treating the vortex sectors w and the coupling configurations J on equal footing is also experimentally motivated. Given sufficient site addressability, the local control of the couplings through external laser fields is how one could perform vortex creation and adiabatic transport in the proposed optical lattice experiments [40,41]. ‡ Formally one should also manipulate the local values of K. However, when K is small compared to the couplings J ij , the magnitude of the three spin term does not significantly affect the physics of the vortices and thus controlling the couplings J ij is sufficient. § Introducing a complex phase is not possible, because the couplings have to be real at all times.

System of interest
To study the physics of the anyonic vortices, we employ a large finite toroidal system of 1200 spins defined by L x = 30 and L y = 20 and consider the system initially in the vortex-free sector by setting u ij = 1 on all links (any other gauge giving rise to w p = 1 on all plaquettes would also do). While this sector of the honeycomb lattice model supports several topological phases where the vortices can behave as either Abelian or non-Abelian anyons [10], we concentrate here only on the latter. It is characterized by Chern number ν = ±1, which implies that the vortices obey Ising anyon statistics. This phase emerges when K = 0 and the couplings violate the inequalities J α > J β +J γ for all all permutations of α, β, γ = x, y, z. Unless otherwise noted, we achieve this by setting J ≡ Jx Jz = Jy Jz = 1 on all links. This amounts to considering the system in the middle of the non-Abelian phase in the sense that the spectral gap is maximized as a function of the couplings J α . For J < 1 the system is closer to the phase boundary located at J = 1 2 . On top of this background we then imprint different vortex configurations by locally changing the couplings as discussed above.
We retain the magnitude K of the three spin term as a free parameter. Its physical range depends on the way it is introduced in the proposed optical lattice experiments [40,41]. If it is introduced by applying a weak Zeeman field [10], only magnitudes of up to K ≈ 10 −4 can be tolerated before the topological phase breaks down [42]. On the other hand, by engineering it directly one can in principle achieve larger values [43].

Vortex-vortex interactions and the fusion degrees of freedom
In this section we show that in the presence of vortices the spectrum exhibits modes whose energy oscillates and converges to zero with the vortex separation. These modes are argued to correspond to the fusion degrees of freedom of the non-Abelian anyonic  vortices. We find the microscopic dependence of the oscillations and show that they are consistent with Majorana fermions being localized at the vortex cores. Finally, we outline the topological low energy spectrum.

Oscillating vortex-vortex interactions
When two vortices are separated linearly as illustrated in Figure 3(a), the spanned configurations are conveniently parametrized by a continuous vortex separation defined as d s = d + s S . Here d denotes the number of links of equal coupling strength between the vortices and 1 ≤ s ≤ S is the step when the instantaneous coupling strength on the (d + 1)th link is 1 − 2s S . At integer values of d s the vortices are pinned exactly to the plaquettes, while for non-integer values, as illustrated in Figure 2(a)-(d), they are interpreted being located at some intermediate position. As d s → 0, the vortices are brought to the same plaquette. We refer to this as fusing the vortices. Figure 3(b) shows the spectral evolution as a function of the vortex separation d s . As observed in [22,44], in the presence of a pair of vortices the energy of the lowest lying mode, E ds ≡ E ds 1 , converges to zero with vortex separation. The finite energy at small separations is interpreted as being due to short-range vortex-vortex interactions. Higher energy modes are insensitive to it and remain at constant energy ∆ f = E ds 2 that corresponds to the spectral gap. The continuous nature of the transport, however, reveals the oscillations in the energy E ds that have been qualitatively predicted in the context of p-wave superconductors [24]. Due to particle-hole symmetry, the oscillating modes come in ±E ds pairs, but the crossings are genuine. In general, the oscillating energy can be written as 1 the coherence length ξ converges to a constant, while λ acquires vortex separation dependence (this appears as an increase due to λ being calculated from fixed nodes at small d s for which the wavelength first increases), before also converging to a J independent constant.
where non-universal constants λ(J, K) > 0 and ξ(J, K) > 0 depend on the couplings and parametrize the frequency of the oscillations and the convergence of the energy, respectively. We are especially interested in the latter as it gives the characteristic coherence length for particular values of the couplings. Only when d s ≫ ξ, one expects the low energy physics to be fully described by the topological Ising anyon theory. The parameters λ and ξ can be extracted from the plot of ln(E ds ) as shown in Figure  4(a). The linear fit with negative slope confirms the exponential convergence of the energy with vortex separation, and distance between successive dips gives the half of the wavelength. By performing similar linear fits for a range of J and K, we obtain Figure 4(b) that shows ξ and λ as functions of the physical parameters. They show very different behavior in the small and large K limits, which roughly occur when K 0.1 or K 0.1, respectively.
Let us consider first the K 0.1 regime, that is experimentally the more relevant one. There the behavior of both the coherence length and the oscillation wavelength is simple. We find ξ decreasing smoothly as K −1 , with the J dependence becoming significant only for K 0.01. This behavior can directly be attributed to the ∆ ∼ K scaling of the spectral gap, which implies that the coherence length scales in general as ξ ∼ ∆ −1 . For instance, when J = 1 the gap is given by ∆ = 6 √ 3K (see Figure 7(b)), which provides an excellent approximation of the K 0.1 region in Figure 4(b). While the coherence length is predominantly controlled by K, the oscillation wavelength λ depends only on J in this regime. This follows from λ being proportional to the inverse Fermi momentum [24]. To be precise, here it is proportional to the difference of the two Fermi point momenta, i.e. 2π λ ∼ |p + F − p − F |. For small K the Fermi momenta are insensitive to the three spin term, while they move away from each other, i.e. |p + F − p − F | increases, when one approaches the phase boundaries (J → 1 2 ) [32]. Thus the observed decrease in λ with decreasing J is consistent with this picture.
The behavior in the K 0.1 regime is dominated by the behavior of the spectral gap. It converges to a J dependent constant value somewhere in the range 0.1 K 0.2 and can not be increased further by increasing K [22]. When this occurs ξ becomes also a constant on the order of the lattice spacing. On the other hand, the oscillation wavelength λ enters an intermediate regime where it becomes first dependent on the vortex separation d s . In Figure 4(b), where the wavelength is calculated from the small d s behavior, this corresponds to the regime where λ first increases and then decreases. When K is increased further, the large d s behavior takes over and λ converges to the value of two lattice constants. This means that for sufficiently large K the oscillations are essentially absent when vortices are pinned to plaquettes (integer values of d s ). This behavior is in agreement with the results for p-wave superconductors, where changes in the oscillations are known to occur both in the vicinity of phase transitions and for suffiently large gaps [26,25]. Even though the gap saturates in the honeycomb lattice model, the topological phase remains robust even for K ∼ J [42]. We interpret the changes in ξ and λ being due to the specific microscopics of the honeycomb lattice model that can modify the non-universal properties of the topological phase. While both the small and large K regimes support the same non-Abelian anyons, they can therefore exhibit very different microscopic signatures depending on the physical parameters. Understanding the microscopics that control these signatures is relevant, for instance, to anyon-anyon interaction driven phase transitions [29], as well as to any topological quantum computing schemes where read-out of information is based on detecting the energy shifts of the topological states [24].
3.1.1. Rotational anisotropy and transversal transport To characterize the fusion mode energy E ds for arbitrary relative vortex locations, we consider two further transport processes. First, let us consider the spectral evolution for the transport shown in Figure  5(a) for various orientations of the lattice. If we assume the vortex at the lower left corner to be located at the origin, we find energy behaving identically in six sectors bounded by the directions ±(n x + n y ), ±(n x − 2n y ) and ±(2n x − 2n y ). Figures 5(b) and 5(c) show that within these sectors one can define parallel directions along which the oscillations have always the same wavelength with the nodes coinciding. The only difference is the overall exponential damping with the vortex separation when the transport is carried out further away from the origin. This sector structure contrasts with the rotational isotropy in p-wave superconductors [45]. The difference is again in the microscopics, which this time derives from to the underlying honeycomb lattice.
Finally, we consider the fusion mode energy when vortices are transported transversally with respect to each other as illustrated in Figure 6(a). As the oscillations   Figure 6. (a) A parametrization for transversal transport. Starting with a configuration of two vortices separated linearly by x plaquettes, we define the transversal direction to be −n x + n y . d x s denotes the distance from the initial configuration. (b) For an even x the energy splitting may oscillate, but one always arrives at the same sign of the splitting. For an odd x the sign always changes such that the energy splitting is exactly zero at the sector boundary.
are known to behave identically along the sector centers, one expects the energy to be reflection symmetric with respect to the boundaries. This is indeed the case as shown in Figure 6(b) for various distances x from the origin, but we also find that the sign of the energy may change. For an even x the energy may oscillate, but one always arrives at the same sign, while for an odd x the sign always changes such that the energy is exactly zero at the sector boundary. This implies that one can arrive at different conclusions on the energy depending on the path chosen to transport the vortices. As all physical vortex states have to be gauge invariant, i.e. not depend on the path [10], this sets limits on simulating vortex transport through manipulating the couplings. Only tuning the couplings along paths without an ambiguity simulate the evolution of the physical states. This can be achieved by considering always the shortest path, which corresponds to tuning only the J y and J z couplings. The restriction to two out of the three link types arises also naturally in the fermionization techniques that do not require additional gauge symmetrization [36,38].

The low energy spectrum and the fusion degrees of freedom
It has been argued in [22] that the modes having the energy behavior (5), to be referred to here as fusion modes, should be identified with the anyonic fusion degrees of freedom. We explicitly verify this by plotting in Figure 7(a) some of the lowest lying states in the two vortex sector with respect to the lowest lying states in the vortex-free sector (w p = 1 on all plaquettes). At large d s there are two exponentially degenerate states with energies E ds 0 and E ds 0 + E ds that differ by the occupation of the fusion mode. Here Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model 12 E ds 0 denotes the energy of the state with unoccupied fusion mode. We define it by as the state with the lowest energy when vortices are nearby. It follows from this definition, that when the vortices are brought closer and finally fused, the energy of the state with the unoccupied fusion mode first oscillates and then evolves smoothly to E 0 , the ground state of the vortex-free sectors. On the other hand, the energy of the state with the the occupied fusion mode becomes ∆ f , the first excited free fermion state on the vortex-free sector. In the language of topological field theories [6,10], the vortex-free ground state, the vortices and the fermions carry a topological quantum numbers 1, σ and ψ, respectively. They satisfy the fusion rule σ × σ = 1 + ψ, which tells that there are two alternatives how a pair of vortices can behave when they are fused: they can either annihilate to the vacuum or leave behind a fermion. Figure 7(a) demonstrates that this fusion degree of freedom is exactly that of the occupation of the fusion modes .
We also observe that when the vortices are fused, the vacuum channel is energetically favoured. This contrasts with the behavior of quasiholes in the Moore-Read state, for which one finds the ψ fusion channel to be energetically favoured [23]. These different results highlight the different microscopics of the models, that can dramatically affect the non-universal behavior of the topological phases.
While the honeycomb lattice model does not accommodate analytic treatment of the fusion modes, their origin can be understood in the context of p-wave superconductors to which the honeycomb lattice model can be mapped [35,36,38]. There one can explicitly show that vortices bind unpaired massless Majorana fermions [9], whose wave functions exhibit exponentially damped oscillating tails [45]. The width of the wave functions is controlled by the fermion gap ∆ f that can be viewed as a height of a potential well confining them to the vortex cores ¶. When two vortices are nearby each other, the overlap of these tails results in a finite tunneling amplitude between vortex cores, which in turn lifts the degeneracy and gives rise to an oscillating finite energy [24]. Our results of the oscillating fusion mode energies are consistent with this picture. While the underlying lattice makes it hard to visualize the oscillations in the wavefunctions, by taking appropriate linear combinations of the fusion modes one can find Majorana modes that have support only on the sites around an isolated vortex [46]. In Section 4 we will employ this picture of localized Majorana modes to construct an effective lattice model for studying degeneracy lifting in many vortex systems.
Strictly speaking there is no ground state degeneracy in the two vortex sector due to fermionic parity conservation. However, the same study could have been carried out with another vortex pair in the system, in which case it is possible to define degenerate states of same fermionic parity. As long as these additional vortices were far away from the other, this would not affect the spectrum and the presented analysis would apply identically.
¶ This is also in agreement with the general theory that says that the fusion channel degeneracy lifting can always be understood as virtual exchange of a fermionic ψ quasiparticle [21]. As ∆ f gives their minimal energy, such processes become less probable as the fermions become more massive.

Energy gaps and stability
When the vortices are well separated, the topological low energy spectrum of Figure  7(a) can be characterized by the following energy gaps: ∆ f and ∆ 2v are the fermion and vortex gaps, respectively, that describe the minimal energy required to excite a free fermion mode and a pair of non-interactiong vortices. V is a local potential that favours the vortices to be located at the plaquettes. To study their physical parameter dependence, we plot them in Figure (7)(b) as functions of the effective magnetic field K. We observe a crossing in the magnitudes of the energy gaps such that for ). This means that while for all values of K the vacuum channel is energetically most favoured, for K < 0.1 the energy of the system can also be lowered by fusing the vortices into fermions. For large effective magnetic field (K > 0.2) the energy gap becomes a constant, the exact value depending only on J [22] (∆ f ≈ 2.0 for J = 1), and can no longer be made larger by increasing K. The periodic "hopping" of the energies as a function of d s has not been observed before. The minima always occur for integer values of d s , i.e. when vortices are pinned to plaquettes, whereas the maxima occur at d S/2 , i.e. when the transported vortex occupies a composite plaquette twice the size of a regular plaquette. We interpret this energy required to move a vortex to an adjacent plaquette as a local potential of magnitude V , defined asymptotically by (9). As shown in Figure 7(b), it varies only slightly with K and can thus be effectively treated as constant for a given uniform coupling configuration J. Therefore, the vortices can be viewed as living on a uniform periodic background potential whose minima coincide with the sites of the dual triangular lattice.
Together the energy gaps ∆ f and ∆ 2v and the potential V give a measure of the stability of the topological low energy theory against local perturbations. When the perturbations are weaker than ∆ f and ∆ 2v , the spontaneous creation of both fermion and vortex excitations is suppressed. When they are also weaker than V , any existing vortices in the system will remain stationary. ∆ f describes the stability of the topological phase itself, whereas ∆ 2v and V describe stability of the vortex sectors within it.

Degeneracy lifting in many vortex systems
In this section we study degeneracy lifting in many vortex systems. Employing an effective lattice model of Majorana fermions, we show that the hybridized spectrum of the fusion modes can be understood in terms of bi-partite interactions between the vortices.

Free Majorana fermions on a lattice
In the previous section we analyzed the spectral evolution in two vortex systems. A simple generalization is to consider what happens when a single vortex is dragged away from a bunch of vortices. When this is performed for the four vortex configuration shown in Figure 8(a), one still finds qualitatively similar behavior. Figure 8(b) shows that one of the fusion modes to oscillates and converges to zero energy, while the other remains at finite energy. In other words, one of the two fusion modes is delocalized while the other remains localized at the three vortices. However, comparing to Figure 3(b), the energy of both fusion modes, including the amplitude and the wavelength of the oscillations, is modified due to the presence of additional vortices. As more vortices are added, the fusion modes will in general hybridize a new energy band within the spectral gap [32]. The degeneracy lifting due to tunneling Majorana fermions suggests that this energy band can be modelled by Majorana fermions tunneling on a lattice defined by the vortex locations [31]. Such a model possesses two properties that we have already observed on the honeycomb lattice model: it has particle-hole symmetry giving rise to double spectrum and it exhibits exact zero energy states for an odd number of sites. The most general Hamiltonian one can write is where t l are tunneling amplitudes between sites of separation l, s l ij are local Z 2 gauge degrees of freedom on link (ij) and γ i are Majorana fermions satisfying γ i † = γ i and {γ i , γ j } = 2δ ij . The Hamiltonian can be defined on arbitrary lattice geometries, which in our case is fixed by the positioning of the vortices. As the interactions in the honeycomb lattice are exponentially suppressed, we simplify the effective model by considering only the three shortest range interactions of range l = 1, √ 3, 2. We refer to them as t 1 , t √ 3 and t 2 interactions, respectively, that couple sites as illustrated in Figures 9(c) and 9(f). The configurations s l = {s l ij } are called gauge degrees of freedom, because the spectrum does not depend on them explicitly. Instead, the model depends on the introduced effective flux, which we define on plaquette (ijk), i, j and k being the three corners of a triangular plaquette enumerated counter-clockwise, by This evaluates always to either ± π 2 . The set of fluxes on all plaquettes is referred to as a flux sector.
The free parameters of the model are the couplings t l and the flux sectors. Our objective is to study how they are to be fixed such that the spectrum of the Majorana model will approximate the spectrum of the fusion modes. To compare quantitatively the spectrumĒ n = (−E n/2 , −E n/2−1 , . . . , E n/2−1 , E n/2 ) of the fusion modes from an nvortex system to the spectrumǭ n = (−ǫ n/2 , −ǫ n/2−1 , . . . , ǫ n/2−1 , ǫ n/2 ) of a corresponding n-site effective model, we introduce the deviation measure For F = 0 the spectra match exactly. To systematically study the many vortex systems, we restrict to considering homogenous vortex configurations, denoted by C d Nx×Ny , consisting of N x vortices in direction n x and N y vortices in direction n y , all separated by at least d links. In particular, we will consider chain (N y = 1) and ladder (N y = 2) configurations that are illustrated in Figures 9(a) Figure  9(c).
Our strategy for studying the many vortex degeneracy lifting is as follows. We set again globally J x = J y = J z = 1, but as before leave the effective magnetic field strength K as a free parameter. Considering then vortex chain and ladder configurations, we find the fusion mode spectraĒ(K) for the range 0 ≤ K ≤ 0.3. Since the fusion mode energy E ds in two vortex systems can be understood as arising due to tunneling between vortex cores, its magnitude at a fixed vortex separation d s = l should correlate with the tunneling amplitudes between sites of separation l. We assume the simplest possible correlation and use the ansatz for the effective model couplings corresponding to C d Nx×Ny vortex configuration. For the flux sectors, on the other hand, there is no obvious a priori correlation with the underlying model. We assume only that topologically equivalent effective models, i.e. ones differing only by the sparsity d, should be defined on the same flux sector. To find the one providing the correct effective description, we need to consider the effective spectrumǭ(K) =ǭ [t l (K)] over all the possible different flux sectors. 4×1 and sparse C 2 4×1 four vortex chains, respectively. In the absence of also t 2 interactions, there are no non-trivial plaquettes and thus the spectrum does not depend on the flux sector. We find that the t 1 interactions alone provide an excellent approximation, which, as quantitatively shown in Figures 10(c) and 10(d), corresponds to at least F < 0.02 and gets better with increasing K. This is in agreement with longer range tunneling being more suppressed for larger K. We note also that in Figure 10 and (b) the sparse C 2 4×1 chain configurations. The red dashed lines denote fusion mode spectraĒ, the black dash dotted lines the bi-partite fusion mode energies E ds corresponding to t 1 interactions, and the squares and circles the effective model spectra with (on the alternating flux sector) and without the t 2 interactions, respectively. The deviations F for (c) the dense chains C 1 N ×1 and (d) the sparse chains C 2 N ×1 for different flux sectors. The circles correspond to t 1 effective model only, while the crosses and squares correspond to uniform and alternating flux configurations, respectively. The approximation by the effective model compared to the pure t 1 model can only be improved when employing the alternating flux sector. The region K < 0.03 is omitted due to finite size effects.
they are either negligible due to t 2 ≪ t 1 (sparse d = 2 chains) or the bi-partite splitting (13) overestimates their strength due to collective effects (tight d = 1 chains). As these observations apply also to longer vortex chains, we conclude that the t 2 interactions are physically relevant in the small K regime of the vortex chain systems. The necessary flux sector is the one having ± π 2 flux alternating on the plaquettes consisting of two t 1 -links and a single t 2 -link.

Vortex ladders
The main difference between the vortex chains and the ladders is the additional presence of the t √ 3 interactions. However, let us start again with and (b) the sparse C 2 2×2 vortex ladder configurations. The red dashed lines denote the fusion mode spectraĒ, the black dash dotted lines the fusion mode energies E ds corresponding to t 1 interactions, and the circles and squares the t 1 effective model spectra on the uniform and the alternating flux sectors, respectively. The deviations F for (c) the dense C 1 N ×2 and (d) the sparse C 2 N ×2 ladders for different flux sectors. The circles and squares denote flux sectors as above. The region K < 0.03 is omitted due to finite size effects.
the nearest neighbour effective model which has now two distinct flux sectors: the uniform sector with Φ 126 = Φ 256 = Φ 235 = Φ 345 or the alternating sector with Φ 126 = Φ 235 = −Φ 256 = −Φ 345 (see Figure 9(f)). Figures 11(a) and 11(b) for the smallest dense and sparse ladders, respectively, show that a qualitatively good fit can be obtained in the uniform flux sector, whereas in the non-uniform sector the spectrum is gapless. This is reflected in the deviation measures, shown in Figures 11(c) and 11(d), where the uniform flux sector clearly provides a better approximation than the non-uniform one. Therefore, we conclude that having a π 2 -flux on all plaquettes consisting of t 1 links only is another necessary condition for the effective Majorana model to approximate the full model.
In general the approximation by the nearest neighbour model is not as good for the vortex ladders as it was for the vortex chains. However, this error is in general (a) (b) (c) Figure 12.
systematic and can be compensated by scaling the tunneling amplitudes as t l = aE dl for some K indepedent constant a. While the systematic study of such dressing effects is beyond the present work, one can imagine them occurring due to quantum interference of many Majorana wave functions. Such effects are also likely to play a role when the t √ 3 and t 2 interactions are switched on. In this case the flux needs to be defined also on plaquettes consisting of two t 1 -links and a single t √ 3 -link, and on plaquettes with one link of each type. For instance, in Figure 9(f) these correspond to the fluxes [Φ 156 , Φ 125 , Φ 245 , Φ 234 , Φ 356 , Φ 362 ] and [Φ 136 , Φ 135 , Φ 462 , Φ 463 ]. We find that the approximation can in general be improved for some K ranges by introducing the longer range interactions, but unlike in the vortex chains we find no single flux sector that would systematically improve it in the physically relevant small K regime. We attribute this to the studied vortex arrays being too small for such genuinely two dimensional effects to become transparent. We leave their systematic study for future work [47].

Effective flux sectors and the underlying π-flux vortices
We discovered that having a π 2 -flux on every plaquette consisting of only t 1 -links is a necessary condition for the effective Majorana model to approximate the full one. Also, in the presence of t 2 interactions, the approximation for vortex chains could be improved only when one imposed an alternating ± π 2 -flux pattern on the plaquettes consisting of two t 1 -links and a single t 2 -link. We postulate that these conditions are connected to the underlying honeycomb lattice model as follows. The anyonic vortices there are π-flux vortices. Let us assume the flux to be uniformly spread across the hexagonal plaquette a vortex occupies and the Majorana mode to be bound exactly at its center. Then, as illustrated in Figure 12(a), a triangular dual lattice plaquette consisting of three t 1 -links will enclose exactly π 2 -flux arising collectively from the three vortices. On the other hand, the plaquettes consisting of a single t 2 link and two t 1 links, as illustrated in Figure 12(b), span no area and should thus have no flux on them. While trivial flux can not directly imposed on the Majorana model, we interpret the alternating ± π 2 -flux sector effectively providing it. We postulate that when writing down a general effective Majorana model for interacting π-flux vortices, the flux on every possible effective model plaquette should be chosen to be consistent with the enclosed vortex flux from the underlying microscopic model. For instance, as illustrated in Figure 12(c), the plaquettes consisting of a single t √ 3 -link and two t 1 -links should have also π 2 -flux. While the concept of enclosed vortex flux is equivalent of the net phase in vortex arrays [31], we emphasize the conceptual simplicity of our prescription in relating the microscopic model and the necessary conditions it imposes on the effective model. The verification of this picture and the development of a thermodynamic Majorana model is a subject of future research [47].

Conclusions
We have studied the interacting non-Abelian Ising anyons emerging from Kitaev's honeycomb lattice model. By simulating continuous vortex transport [30], we uncovered the characteristic oscillations in the bi-partite degeneracy lifting and characterized their physical parameter dependence. While the oscillations had been discovered earlier for both p-wave superconductors [24] and the Moore-Read state [23], these calculations involve mean fields and trial wave functions, respectively. Our results demonstrate the oscillations in an actual microscopic model. Due to algebraic similarity, we expect them to be present also in the variations of the honeycomb lattice model [11,12,13,14,15]. Furthermore, we extended the results of [22] by demonstrating unambiguously the fusion degrees of freedom and by obtaining the energy gaps characterizing the stability of the topological low energy theory.
In the second part we studied degeneracy lifting in interacting many vortex systems. Employing the picture of localized Majorana modes, we wrote down a Majorana fermion model on a finite lattice whose sites coincide with the vortex locations. By comparing its spectrum to that of the fusion modes of a corresponding vortex configuration, we were able to find necessary conditions for the spectra to match. This amounted to finding relevant flux sectors of the effective model, which we interpreted as arising from the πflux vortices of the underlying full model. These results are consistent with and extend the earlier mean-field studies in the context of p-wave superconductors [31]. Further, we found that the energy splitting due to vortex-vortex interactions could be employed to high accuracy as the effective tunneling amplitude of the Majorana fermions. This confirms that bi-partite tunneling is predominantly responsible for the degeneracy lifting also in many vortex systems. If many body effects are present, they can be included as dressed or longer range couplings. Our results pave the way for constructing an effective model for interacting anyonic π-flux vortices at the thermodynamic limit [47]. Such model can be employed to understand phase transitions due to vortex lattices, [32], and test the general theory of anyon-anyon interaction driven phase transitions in the context of the microscopic honeycomb lattice model [29,47].