Majorana fermions on a disordered triangular lattice

Vortices of several condensed matter systems are predicted to have zero-energy core excitations which are Majorana fermions. These exotic quasi-particles are neutral, massless, and expected to have non-Abelian statistics. Furthermore, they make the ground state of the system highly degenerate. For a large density of vortices, an Abrikosov lattice is formed, and tunneling of Majorana fermions between vortices removes the energy degeneracy. In particular the spectrum of Majorana fermions in a triangular lattice is gapped, and the Hamiltonian which describes such a system is antisymmetric under time-reversal. We consider Majorana fermions on a disordered triangular lattice. We find that even for very weak disorder in the location of the vortices localized sub-gap modes appear. As the disorder becomes strong, a percolation phase transition takes place, and the gap is fully closed by extended states. The mechanism that underlies these phenomena is domain walls between two time-reversed phases, which are created by flipping the sign of the tunneling matrix elements. The density of states in the disordered lattice seems to diverge at zero energy.


Introduction
Exotic states of matter are among the most intriguing topics in the field of condensed matter physics. One class of these exotic states are the two-dimensional (2D) systems in which quasi-particles follow non-Abelian quantum statistics [1]. The search for such systems is driven both by their unique properties and by their potential application to topological quantum computing [2].
In a non-Abelian state the ground state is degenerate when quasi-particles are present, and the degeneracy increases exponentially with the number of quasi-particles. Perhaps the simplest way of obtaining such a degeneracy is by having quasi-particles that carry Majorana fermionic excitations [3]. Majorana fermions (MFs) are expected to appear as zero-energy excitations in the cores of vortices in a layered p x +ip y superconductor -such as proposed for Sr 2 RuO 4 [4], in 2D systems that can be mapped onto such a superconductor -the ν = 5/2 fractional quantum Hall state [5], on the surface of a topological insulator that is in proximity to an s-wave superconductor and an insulating ferromagnet [6], and at hybrid structures of semiconductors and superconductors [7]. Furthermore, they are expected to form in 1D systems in proximity to s-wave superconductors [8], and in several other systems [9,10].
A variety of experiments have been proposed in order to probe the predicted MFs, based on the their unique properties. To mention a few examples: zero-energy excitations can be observed by performing STM measurements at the vortex core [11,12]; the degeneracy of the ground state affects the thermodynamical properties [13,14]; non-locality of an electron in a Majorana state has a signature in tunneling between two vortices [15,16]; and interferometric experiments are able to probe the non-Abelian statistics [2,17,18,19].
The suggested thermodynamics and interferometry measurements are based on the unique many-body properties of the MFs, but assume that the MFs are localized at theirs positions. The MFs appear at vortex cores, which are expected to rearrange as an Abrikosov lattice at high enough density. The lattice order and the small tunneling amplitude between neighboring vortices remove the ground state degeneracy and form a band of low-energy excitations. On one hand, this band may serve to probe the existence of the MFs. On the other hand, it may conceal signals of suggested manybody measurements, especially controlled adiabatic processes, such as interferometry.
Some previous works have analyzed clean periodic square, triangular [20] and honeycomb [21] lattices of MFs, and found the electronic conductivity associated with tunneling between vortex cores. Other works have considered some aspects of random square [10] and honeycomb [22] lattices. In this paper, we consider disordered triangular lattices. This lattice breaks time-reversal symmetry: it may be mapped onto a tightbinding model of electrons on the same lattice with each plaquette being pierced by ±1/4 magnetic flux quantum. The sign of the flux is determined by the sign of the tight-binding coupling term. The spectrum of the perfect lattice is gapped.
Our main finding is that disorder in the sign of the tunneling amplitude between neighboring sites creates a peak in the density of states (DOS) at zero energy. We show that at weak disorder zero-energy localized sub-gap states are created, and coupling between these states creates a DOS close to zero energy. At strong disorder the sign flipping creates domain walls between regions of the two time-reversed phases, and chiral modes appear along these walls. When the signs are random, we find the domains to be narrow and the spectrum of excitations to be characterized by strong dependence on the geometry of the walls. Plausibly, that is the source of the zero-energy DOS peak in the highly disordered system. The transition between these two limits is a percolation phase transition, taking place when the probability of flipping a sign is around 0.15. We find the DOS to diverge at zero energy at intermediate and strong disorder. The paper is organized as follows: In section 2, we define the MFs lattice model, and show the numerical DOS for a disordered lattice. In section 3, we examine how low-energy excitations may emerge from localized modes in the limit of weak disorder and from extended interfaces in the limit of large p. Section 4 discusses the dependence of the excitation spectrum on the strength of the disorder. Section 5 summarizes the results and compares them to previous works.

A disordered triangular lattice
Our interest here is in MFs that form a triangular lattice. An MF is defined as a selfadjoint operator γ = γ † that satisfies fermionic anti-commutation relation {γ i , γ j } = δ ij . A standard complex fermionic operator can be constructed from MFs by superposing an even number of them. For example ψ n = (γ i + ıγ j )/ √ 2 for two MFs satisfies the standard anti-commutation relations {ψ n , ψ † m } = δ nm and {ψ n , ψ m } = 0. We consider MFs that are solutions of the Bogoliubov de-Gennes (BdG) equation in the presence of vortices in the superconductor. Each vortex in the superconductor carries a signle MF, whose wavefunction is exponentially localized around the vortex's core. Overlaps of the MF wavefunctions of neighboring vortices result in tunneling matrix elements. Vortices in superconductors and quasi-particles in clean quantum Hall systems are arranged on a lattice; thus the MFs are also arranged as a lattice, and their dynamics may naturally be described by the tight-binding model.
The particle-hole symmetry of the BdG equation implies that a single MF is a zeroenergy solution. Hence, in the tight-binding Hamiltonian the on-site energy of the MFs is zero, and there are only hopping terms t ij . Moreover, for the Hamiltonian to be Hermitian, t ij γ i γ j = (t ij γ i γ j ) † = −t * ij γ i γ j , implying a purely imaginary hopping term. By assuming discrete symmetry of the lattice to translations, we can write the simple Hamiltonian where ij are nearest neighbors and s ij = −s ji = ±1. Any element s ij is gauge dependent, because the γ i operators are defined up to an overall sign. However, the product of s ij 's along a path creating a closed loop is gauge independent. It has been shown in [20] that for any lattice whose plaquette is a polygon of n vertices, the product of s ij around each plaquette corresponds to the plaquette enclosing n/4 − 1/2 flux quanta. Therefore the product of the hopping terms along the bonds that create the plaquette is −ı n t n , and the product of the s ij 's along this path is −1 (the direction of the path is chosen to be aligned with the chirality of the order parameter).  (2). An arrow from site i to site j means s ij = 1. Note that the lattice is split into two sublattices A and B, and that the counterclockwise product of the s ij 's around each triangle equals to −1. (b) The corresponding gapped spectrum, with E gap = 1.73.
In particular, in the triangular lattice each plaquette encloses a quarter of flux quantum, and the product of the hopping terms equals ıt 3 , revealing a time-reversal anti-symmetry. Moreover, the unit cell which encloses a flux quantum is a parallelogram of four neighboring triangles, leading to the lattice being composed of two sublattices. One possible gauge is illustrated in figure 1(a), which also depicts the two sublattices A and B. In this gauge Let us assume a clean periodic lattice of L 1 × L 2 sites, with L 2 even. We can define the Fourier transformΓ a,k = i e −ıkx i γ a,i , where a = A, B and k = 2π(m 1 /L 1 , m 2 /2L 2 ) for m i = 0, ..., L i − 1 (i = 1, 2) [20]. These transformed operators are complex fermions, which satisfy {Γ a,k ,Γ † a,k } = δ kk . Note, however, that for k = (0, 0), (π, 0), (0, π/2) and (π, π/2) they are MFs. By denoting spinorΓ † k = (Γ † A,k ,Γ † B,k ) the Hamiltonian can be expressed as where the Pauli matrices act on the sublattice space. The resulting spectrum has an energy gap of 2tE gap , where E gap = 1.73, as depicted in figure 1(b).
The hopping elements between neighboring MFs are very sensitive to the inter-vortex separation. For example it was shown that for a p x +ip y superconductor [23]: where k F is the Fermi momentum, r is the inter-vortex distance, and ξ is the coherence length, which is usually larger than k F . The exponential decay and the oscillations were found to occur also for the ν = 5/2 case [24], and are likely to appear in all realizations. Therefore small deformations of the Abrikosov lattice of order k −1 F will produce fluctuations in both the amplitude and sign of the hopping terms t ij . In systems with hopping terms that are random both in sign and in magnitude (blue) and in systems where only the signs are random (green) there is a sharp peak at zero energy, while random amplitudes (red) only smear the spectrum of the clean system, but do not close the gap. The DOS is averaged over 10,000 realizations.
We find numerically that random hopping terms t ij = ts ij , where s ij is uniformly distributed in the interval [−1, 1], close the energy gap of the uniform system. The DOS, which appears in figure 2, shows that not only the gap is closed, but a peak emerges at zero energy. We can distinguish between random amplitudes and random signs of the hopping terms, i.e. |s ij | ∼ U [0, 1] or s ij = ±1 in equal probability, respectively. The DOS of these two cases, which are also depicted in figure 2, clearly shows that the zeroenergy peak in the DOS appears only due to random signs, while random amplitudes merely smear the spectrum without closing the energy gap. We are mostly interested in the zero-energy peak of the DOS, and therefore in the following we will focus on the case where disorder appears in the sign of the hopping terms.
The sharp peak that we found numerically for the density of states of finite systems close to zero-energy naturally raises the question of the way the DOS scales with the size of the system in the thermodynamic limit. Increasing the system size L obviously increases the DOS. If the lowest energies decay faster than 1/L 2 , then the zero-energy DOS will increase faster than L 2 , and the zero-energy DOS per unit area would diverge in the thermodynamic limit. Figure 3(a) exhibits the energy of the three lowest states E 1 , E 2 and E 3 , averaged over disorder, as a function of the system size. In the given range a fit to E n ∝ 1/L α yields α ≈ 2.18 for n = 1, 2 and 3, which indicates either a weak power-law or logarithmic divergence of the DOS per unit area.  It is instructive to compare the DOS we find for the case of random nearest-neighbor hopping to that we find for the case of a random Hermitian matrix of imaginary terms. The DOS of the latter is a version of a semicircle distribution, with a Delta peak at zero energy, due to the symmetry of the spectrum with respect to reflection about zero energy [25]. When we add second nearest-neighbor hopping with random signs to the Hamiltonian, the DOS indeed gets closer to the semicircle. Figure 3(b) shows the spectrum for several ratios of t 2 , the amplitude of the next-nearest-neighbor terms, to t. We note, however, that we numerically find the gap to close even for the clean system when t 2 /t ≈ 0.58.
Having established that the DOS of the clean lattice is characterized by an energy gap limited by two singular peaks, and the lattice where the sign of the hopping terms is random is characterized by a zero energy peak above a gapless background, we examine the evolution of the DOS with p, the probability of flipping the sign of a hopping term, which increases from 0 (the clean lattice) to 0.5 (the random sign limit). Figure  4(a) shows how increasing p makes the singular peaks of the clean system spread, and gradually creates sub-gap states.
In the next section we study the way the gap is closed, and especially the appearance of the zero-energy peak. We address these questions by an examination of two limits. In the weak disorder limit, we present the minimal disorder configuration that results in a zero-energy state. In the strong disorder limit, we examine the low-energy states that are formed along domain walls between regions in which the signs of flux piercing the plaquettes are opposite. Figure 4. The disorder-averaged DOS D(E) of a periodic 70 × 70 lattice as a function of the probability p for flipping the sign of a hopping term. For p = 0 the DOS is that of figure 2(a). Two small peaks at E ≈ ±1, that grow linearly with p, belong to localized states which are created around a single isolated flipped hopping term. The zero-energy peak grows slowly, but survive at the disordered system. The DOS is averaged over 100 realizations.

The origin of the zero-energy peak in the density of states
We have seen that the clean triangular lattice of MFs is anti-symmetric under time reversal, and has an energy gap. Depending on the sign of the hopping tunneling elements, the MFs may be mapped onto electrons in a triangular lattice pierced by 1/4 or −1/4 flux quanta per plaquette. These two cases are characterized by two opposite Chern numbers: by writing the k th component of equation (3) as H k = 2th k · σ, it is easy to see that the Majorana band is characterized by a non-trivial Chern number The disorder-induced reversal of signs of hopping terms may revert the sign of the flux in some of the plaquettes in the lattice. The systems would then have islands of one phase separated by lines of interface from the bulk of the other phase. An infinite interface between two phases of different Chern numbers is accompanied by gapless modes [26,27,28]. For finite islands, one may naively expect finite-size quantization to induce a gap in the energy spectrum, and thus low-energy excitations to require large islands. This expectation is only partially valid. As we explain below, for the Majorana lattice, large islands are associated with low-energy excitations at their edges, but small islands may carry low-energy and zero-energy modes as well.
When the probability p of flipping the sign of a hopping term is very small the flipped terms are dilute, and most of them are isolated. A two-triangles island, which is created by flipping the sign of a single hopping term, results in two localized sub-gap states with energies of approximately ±1.1. Theses states are the source of the peaks in the DOS at ±E ≈ ±1.1, and indeed in the limit of small probability the DOS associated with these peaks D(E ≈ ±1.1) ∝ p. Such islands do not, however, contribute to the DOS close to zero energy.
As we show in Appendix A, flipping three hopping terms with a common vertex creates two localized states with energy that is either zero or exponentially small with the system size. We also show that this is the minimal way of creating zero modes. Thus, for small p the zero-energy DOS is dominated by the probability of creating such configurations, and should scale like p 3 . The two modes are, however, split when there is a single flipped hopping term at a distance r from that vertex. The exact splitting depends on the orientations of the bonds with flipped signs, but we found numerically that it can be well approximated by ±E t (r) ≈ ±6e −2r . For E t to be much smaller than the typical finite-size energy splitting, that scales as L −2 , we need r log L. A given radius r encloses approximately 3πr 2 hopping terms around the vertex, and the probability that all these terms are unflipped is (1 − p) 3πr 2 . Thus, the p 3 dependence of the zero-energy DOS is limited to small values of p log −2 L. For example, for a 70 × 70 lattice we get p < 0.01.
As p gets large, the plaquettes in which the flux is reversed connect to one another, and the system is split into domains with opposite fluxes. To understand the excitations spectrum in this limit, we first examine the spectrum of excitations associated with the various possible interfaces of two large domains of opposite fluxes. Then, we examine the statistical distribution of such interfaces as a function of the flipping probability p.
At the interface between a macroscopic region of Chern number ν = ±1 and the vacuum, a chiral gapless mode must appear [26,27,28]. In order to explicitly find this mode, we first replace η 1 and η 2 in equation (2) byx andŷ for simplicity, and denote the lattice sites by x. We assume rotational symmetry alongŷ, and define the operators Γ a,kx = y e −ıky γ a,x with a = A, B and k = πm/L y , where m = 0, ..., L y − 1. Each Γ a,kx represents a wavefunction that is localized atx but extended atŷ. With these operators the Hamiltonian becomes Given a right edge at x = 0, the edge eigenmode of (6) for For k = 0 it reduces to (1, 1)Γ 0,0 . The dispersion of the edge mode is E k = 2 sin k, with E k ≈ 2k = 2πm/L y at low energy. Note that the wavefunction is localized in thex direction with a localization length of the order of 1/ log(L y /πm) ∼ 1/ log E, which depends weakly on the energy and the system size.
A triangular lattice with one periodic coordinate can have either flat or zigzag edges, depending on the way the periodic boundary condition is taken, as illustrated in figure  5. The Hamiltonian (6) describes a flat edge. In a zigzag edge, we also get a chiral edge mode, with the dispersion E k ≈ √ 2k and b k ≈ (3 − 2 √ 2) − (8 − 11/ √ 2)k 2 , which gives a localization length that is approximately independent of energy, at low energies. We can see that in spite of the differences, both edges show a similar behavior of linear dispersion and exponentially localized wavefunction. Henceforth we will assume the periodicity of equation (6). An interface, or domain wall, between the two phases of opposite Chern numbers is expected to hold two chiral edge modes, similar to the single mode that exists along the interface between any of the two phases and the vacuum. We denote the unflipped phase by P and the flipped phase by N, representing their positive and negative fluxes respectively. A periodic lattice with an interface between a P-phase and an N-phase may have two kinds of domain walls that preserve the translational symmetry: flat and sawtooth, as depicted in figure 5.
A flat domain wall is created when σ y → −σ y in the hopping terms for x > 0. The two chiral modes of such a domain wall have the same dispersion and amplitudes b k as the flat edge mode, where in the N-phase the amplitudes are b * k , as expected from the time reversal. Only now one mode has non-vanishing amplitudes at x = 2n, while the other mode has non-vanishing amplitudes at x = 2n + 1. A sawtooth domain wall is a result of flipping the σ x 's. It shows essentially a similar behavior to the flat domain wall, but with E k ≈ 1.95(k ± k z ), where k z = arccos( √ 5/2 − 1/2). Again, we see that the geometry of the domain walls has only a small effect on the low-energy physics.
If an 'island' of a roughly circular shape of the N-phase is surrounded by the P-phase lattice, then the two chiral modes will be localized along the perimeter of the island, and their dispersion will be E m ∼ m/M , where M is the perimeter of the island. Such an island will close the energy gap only if it becomes macroscopic. In a disordered system where the sign of the hopping term s ij is random, islands of the N-phase will be created in various shapes and sizes. A way to parameterize the morphology of an island is by the perimeter-area relation. In approximately rounded shapes M ∝ A 1/2 , where A is the area of the island, while in a narrow strip M ∝ A. Figure 6(a) depicts the distribution of the perimeters as a function of the areas of random islands in a 30 × 30 lattice, where the probability of flipping a sign is 1/2. The area of an island is the number of N-triangles that compose the island, where two N-triangles are defined to belong to the same island if they share at least a vertex. The perimeter is composed of the P-triangles which have a common vertex with the island. For islands with A > 10 the mean perimeter can be excellently approximated by M (A) ≈ 2.34A + 15.6, which means that the islands look like snakes or gossamer. The smallest island that may be created by flipping the sign of a single hopping term is composed of two triangles, and has a rhombus shape. The simplicity of the shapes is not preserved for larger islands. Flipping two neighboring hopping terms results in two flipped N-triangles, but those share only a vertex. Larger islands are typically of complicated shapes, composed of rhombi and triangles, as illustrated in figure 6(b). Complicated islands provide complicated dispersion and DOS, as we will demonstrate. In order to get a feeling to the resulting dispersions, we now focus on three types of narrow islands that preserve the translational symmetry, i.e. strips of N-triangles.
The method we use for approximating the low-energy dispersion of the strips is first to find zero-energy modes with k = 0, and then extract the dominant k-dependence employing first order perturbation theory. For the clean system the Hamiltonian (6) may be recast as For any strip s, this Hamiltonian has to be modified to the Hamiltonian H s , which may be decomposed similarly to H s 0 and δH s k . The k = 0 zero modes of H s will be denoted by Γ s i = Γ s † i . They satisfy [H s 0 , Γ s i ] = 0. In the cases we will explicitly consider here, there are two zero modes per strip. The matrix elements of δH k between these zero modes give the leading order of the k-dependence. Therefore after finding Γ s I and Γ s II , } The first and simplest strip is the flat strip, which is depicted in figure 7(a), and is created by flipping the sign of the σ z 's in equation (7b) between x = 0 and x = 1. Now H flat 0 has two zero modes: The matrix elements of δH k in the zero modes space are with τ being the vector of Pauli matrices. The resulting corrections to the zero energies are δE flat k ≈ ± √ 5tk. These modes are essentially the domain wall modes, which are modified due to their hybridization, and similarly their DOS is a constant.
The second example is a 'zipper' strip, shown in figure 7(b), which is a result of replacing σ y of Equations (7b) and (7c) by ıσ x for x = −1 and by −ıσ x for x = 0. Two zero modes again appear The third example is the sawtooth strip, which is created by multiplying the hopping terms between x = −1 and x = 0 by σ z , and is shown in figure 7(c). The zero modes are now Unfortunately, [δH k ] saw = −t(1+τ z ) sin k, which means that the first order perturbation theory gives the correction only to the first mode δE saw k ≈ −2tk, while the second mode requires higher orders. An explicit solution of this mode gives δE saw k ≈ k 3 ; thus this soft mode makes the DOS diverge as E −2/3 . We can conclude from this section that along edges and domain walls chiral modes appear with linear low-energy dispersion. These modes close the energy gap only for macroscopic domain walls, which will give a constant DOS. However, when the domains become narrow strips, the dispersions of the chiral modes strongly depend on the geometry of the strip, and their DOS tends to diverge at zero energy. And since most domains which are created by random hopping terms are composed of narrow strips, this may explain the zero-energy peak raising above the uniform background in the DOS.

Percolation phase transition
In the previous section we have seen that chiral modes which are localized along stripeshaped domains may lead to a zero-energy peak in the DOS when the signs of the hopping terms are random. Such modes provide low-energy states (states whose energy scales inversely with the size of the system) only if the length of the domain is of the order of the system size. For hopping terms whose signs are random (p = 0.5) there are such extended domains, since triangles of both phases are distributed all over the lattice. In contrast, for small p only small isolated domains are created.
In this section we address the dependence of the DOS on p. We start by examining the dependence of the size of the islands on p. It is natural to expect a percolation phase transition, where below a critical probability p c only microscopic islands can be found in the lattice. The probability distribution to find an island with area A is then exponential with A, where the typical area depends on p and is independent of the system size L. Above p c a macroscopic island forms; thus the probability distribution is concentrated around the macroscopic value, which of course scales as L 2 .
This expectation is borne out by our numerical analysis. Figure 8 depicts the distribution of the islands area A, normalized to the system's size L 2 , as a function of p for an L = 100 lattice. For large p the distribution appears to be concentrated around a mean value, which was numerically verified to scale as L 2 . In contrast, for small p the distribution is approximately exponential in a manner that was found to be independent of L. The transition takes place at p c ≈ 0.15. The width of the transition between the microscopic and macroscopic distributions, which takes place here in 0.10 < p < 0.15, gets sharper towards p c while increasing L.
Since flipping a hopping term flips the fluxes of two triangles, islands with even values of A are more common than islands with odd values of A. Figure 8 shows that in spite of the different distribution of the even and odd values of A, both are exponentially distributed below p c .  Figure 8(d) depicts, for every p and normalized area A/L 2 , the probability that there are islands larger in area than A/L 2 . The curve for which this probability is 0.5 is the typical size of the largest island for a particular disorder realization. We denote this curve by A max (p). If we measure the flux in the system with reference to the flux of the clean system, then since a triangle with N flux is formed by flipping either one or all three hopping terms forming the triangle, the mean flux threading the lattice is Above p c the macroscopic island, if it exists, is supposed to capture almost all the N flux. Indeed, A max (p > p c ) is excellently approximated by φ .
Thus, above p c almost all the triangles with flipped flux are connected. Figure 9 depicts two low-energy disorder-averaged quantities for various lattice sizes L. The first is E 1 (p) , the energy of the lowest excitation. The second is the zeroenergy DOS D 0 , which we define as the number of states within the energy window 0 ≤ E < E sat , divided by E sat , where E sat = E 1 (p → 0.5) is the disorder averaged energy of the lowest excitation for the p = 0.5 case, extracted from figure 3. In order to compare systems with different sizes, we normalize E 1 (p) by E sat (L), and D 0 by D sat = D 0 (p → 0.5) .
For small p we find E 1 to approach the value of the energy gap, as it should. We find D 0 to be polynomial in p for small p, in accordance with our estimate for the probability of finding isolated regions of zero-energy states. Above the percolation phase transition the mean minimal energy E 1 and the zero-energy DOS D 0 are expected to saturate to a constant value as a function of p. We do indeed observe a saturation of E 1 and D 0 above p c . Surprisingly, figure 9(a) shows that the approach to the saturated value is not monotonic. The minimal E 1 appears at p e ≈ 0.11 < p c , and we find E 1 (p e ) ∝ L −2.88 , as shown in figure 9(c). Moreover, the maximum D 0 also appear at p e . The divergence of the DOS at zero energy seems therefore to be strongest before the percolation threshold.
If the percolation phase transition captures the essence of the low-energy physics, we naively expect the low-energy states to be extended above p c , while those below p c , including p e , to be localized. For any eigenstate of the Hamiltonian Γ E = i a i γ i , the localization properties of the state can be characterized by the participation ratio (PR), defined by P R = ( i |a i | 4 ) −1 . The PR of a perfect metal scales linearly with the system size L 2 , while the PR of an insulator is independent of L. Figure 9(d) depicts the disorder-averaged PR of the lowest excitation as a function of the system size, for various values of p. As expected, above p c the PR scales faster than L, which indicates an extended state, while well below p c the PR is independent on L. However, at p e the PR scales as L. This scaling seems to imply that there is a range of p in which the low-energy states are string-like in the sense that they are macroscopic only along one dimension. Recalling the observed diverging DOS of the edge states along the narrow strips, the divergence of D 0 at p e is understood, even if not the precise exponent.

Summary
In the previous sections we have seen that the zero-energy peak of the DOS for random signs in the nearest neighbor hopping terms may have different sources in two limits of the probability p for flipping the sign of individual hopping terms. For small p the peak is the consequence of the existence of localized zero-energy states and the tunneling between them and other localized states. Therefore it appears on a background of the zero DOS of the energy gap, and it is proportional to p 3 . In the opposite limit of p c ≤ p ≤ 0.5, where p c ≈ 0.15, macroscopic islands of the time-reversed phase cross the  Figure 9. (a) The normalized mean minimal energy E 1 /E sat as a function of the probability p, where E sat = E 1 (p → 0.5) . The percolation phase transition is manifested as the saturation to constant values for p > p c , but surprisingly the lowest energies appear around p e < p c , which means that they belong to localized states. (b) The mean zero-energy DOS D 0 = D(0 ≤ E < E sat ) as a function of p, normalized by its value at p = 0.5, D sat , of the same lattices. The saturation for p > p c is expected, but the DOS gets its maximal values at p e < p c , and decays towards zero only in a polynomial manner. (c) The mean minimal energy at p e as a function of the system size is very well fitted to L −2.88 . (d) The disorder-averaged participation ratio P R of the minimal energy state as a function of the system size L, for various values of p. For p > p c the state is extended, since P R scales faster than L (denoted by dashed line). For p p c the states are localized, since P R is independent of L. At p ≈ p e , P R ∼ L, which implies a 1D string-like state.
lattice, and the gap is closed due to the low-energy modes along the interfaces of these islands with the original phase. Since the dispersion of some of these modes is nonlinear, they are expected to create zero-energy peak in the DOS. Our numerical results indicate a weak divergence of the DOS at zero energy for p = 0.5, and an even stronger divergence somewhat below the percolation transition. In this region the islands seem to be of the form of strings, and the mentioned nonlinear dispersion is likely to be the source the divergence.
It was mentioned above that the Hamiltonian of our system is purely imaginary and anti-symmetric both to charge conjugation and time-reversal. According to the standard classification of disordered Hamiltonians [29] it belongs to class D [30,31]. This class is known to have three distinct phases: thermal metal, thermal insulator and thermal quantum Hall (TQH) insulator [32,33]. In the metal phase, the DOS at zero energy diverges logarithmically with the system size [32]. Several numerical studies of the properties of this class have been carried out using the Cho-Fisher network model [34,35,36], which exhibits all of the phase diagram. It was found that in the TQH phase there is a region below the insulator-metal transition where the DOS at zero energy diverges, but with a nonuniversal exponent, due to the effects of rare configurations of disorder (Griffiths phase) [36]. Recently a model of MFs on a square lattice created in a p-wave superconductor by a random potential has been proposed [10], and the divergence in the metal phase has also been observed numerically there.
Our model is a new realization to the TQH-thermal metal transition, which may be physically realizable. We observe the divergence of the DOS both in the metal regime and in the analogue to the Griffiths regime, and the microscopic understanding of our system provides insights into the physics of these phenomena.
Note added: After completing the preparation of this manuscript, we became aware of a study of a similar system carried out by C Lohmann, AWW Ludwig and S Trebst [37].
swapping labels of two sites, also flips the sign of Pf(S), since it can be performed by B = I (N −2)×(N −2) × σ x , where σ x is the Pauli matrix which acts in the space of the two swapped sites.
Systematic gauge transformations reveal equivalence relations between apparently different lattices. The simplest example is the gauge transformation that swaps between the two sublattices, which is depicted in figure A1(a). In a periodic lattice the number of rows is even (which means that the two sublattices are indeed equivalent), an even number of sites are gauged, and the Pfaffian is unchanged. Rotations by 60 • clockwise and counterclockwise can be also produced by gauge transformations, as depicted in figures A1(b) and A1(c). In principle, such rotations may multiply the Pfaffian by a factor of −1. However, three 60 • rotations result in swapped sublattices, which does not change the Pfaffian, and hence one rotation must leave the Pfaffian unchanged as well. Note that in finite lattices the 60 • rotations may have O(e −L ) corrections, due to change of boundary conditions. All these gauge-dependent properties of the Pfaffian show that it is an unphysical quantity. However, since Pf(A) 2 = det(A), the Hamiltonian H has zero energy states iff Pf(S) = 0. Therefore the way to prove that flipping three neighboring arrows gives zero-energy states, would be to prove that the Pfaffian vanishes for such a configuration.
Let us denote an arbitrary lattice site by 1, and its six neighbors by a, b, .., f , as depicted in figure A1 The physical meaning of Pf(Ŝ 1i ) is omitting the sites 1 and i. We look for relations between lattices with vacancies of two neighboring sites. These relations will allow us to express five of the terms Pf(Ŝ 1i ) in (A.4) in terms of the sixth, say Pf(Ŝ 1a ). Due to the horizontal symmetry of our gauge, it is apparent that |Pf(Ŝ 1b )| = |Pf(Ŝ 1a )|, but since the labeling of the sites is different, they may differ in signs. If, however, we 'push' the a th row and column ofŜ 1b a − b − 1 rows and columns upwards, we get exactlŷ S 1,a . This process involves a − b − 1 labels swapping; thus Pf(Ŝ 1b ) = (−1) a−b−1 Pf(Ŝ 1a ). S 1c andŜ 1d are related by swapping the two sublattices; thus |Pf(Ŝ 1d )| = |Pf(Ŝ 1c )|. By 'pushing' the row and columns, we get again Pf(Ŝ 1d ) = (−1) c−d−1 Pf(Ŝ 1c ). In a similar manner we have Pf(Ŝ 1f ) = (−1) e−f Pf(Ŝ 1e ). Furthermore,Ŝ 1c (Ŝ 1e ) is related toŜ 1a by the anticlockwise (clockwise) rotation, hence |Pf(Ŝ 1e )| ≈ |Pf(Ŝ 1c )| ≈ |Pf(Ŝ 1a )|. By counting the number of sites been gauges, we get Pf( In the clean lattice s 1a = −s 1b = s 1c = −s 1d = −s 1e = −s 1f , as can be seen in figure  A1(d). Thus the terms sum-up, and |Pf(S)| = 6|Pf(Ŝ 1a )|. On the other hand, flipping three of the s 1i 's, gives Pf(S) = 0, i.e. zero-energy modes. Due to the locality of this manipulation, the zero modes are localized around site 1. On the other hand, flipping only one or two terms will not result in zero modes.