Counting Majorana bound states using complex momenta

Recently, the connection between Majorana fermions bound to defects in arbitrary dimensions, and complex momentum roots of the vanishing determinant of the corresponding bulk Bogoliubov-de Gennes (BdG) Hamiltonian, has been established (EPL, 2015, $\textbf{110}$, 67005). Based on this understanding, a formula has been proposed to count the number ($n$) of the zero energy Majorana bound states, which is related to the topological phase of the system. In this paper, we provide a proof of the counting formula and we apply this formula to a variety of 1d and 2d models belonging to the classes BDI, DIII and D. We show that we can successfully chart out the topological phase diagrams. Studying these examples also enables us to explicitly observe the correspondence between these complex momentum solutions in the Fourier space, and the localized Majorana fermion wavefunctions in the position space. Finally, we corroborate the fact that for systems with a chiral symmetry, these solutions are the so-called"exceptional points", where two or more eigenvalues of the complexified Hamiltonian coalesce.


I. INTRODUCTION
Topological superconductors 1 are systems which can provide the condensed matter version of Majorana fermions, as they can host topologically protected zero energy states at a defect or edge, for which the creation operator (γ † E=0 ) is equivalent to the annihilation operator (γ E=0 ). These gapless states obey non-Abelian braiding statistics 2,3 , which can find potential applications in designing fault-tolerant topological quantum computers 2,4 . Although Majorana fermion bound states have not yet been conclusively found in nature, they have been theoretically shown to exist in low dimensional spinless pwave superconducting systems 2,5 , as well as other systems involving various heterostructures with proximityinduced superconductivity which are topologically similar to them [6][7][8][9][10][11] .
Non-interacting Hamiltonians for gapped topological insulators and topological superconductors, in arbitrary spatial dimensions, can be classified into ten topological symmetry classes 12,13 , characterized by certain topological invariants. Here we focus on 1d and 2d Bogoliubov de Gennes (BdG) Hamiltonians with the particle-hole symmetry (PHS) operator squaring to +1, which can be categorized 14 into three classes: BDI, DIII and D.
In our earlier work 15 , we have explored the connection between the complex momentum solutions of the determinant of a bulk BdG Hamiltonian (H BdG ) in arbitrary dimensions, and the Majorana fermion wavefunctions in the position space associated with a defect or edge. We have found that the imaginary parts of these momenta are related to the exponential decay of the wavefunctions, localized at the defects, and hence their sign-change at a topological phase transition point signals the appearance or disappearance of Majorana zero mode(s). Based on this understanding, we have proposed a formula to count the number (n) of the zero energy Majorana bound states, which is related to the topological phase of the system. This formula serves as an alternative to the familiar Z and Z 2 topological invariants 12,13,16 and other counting schemes 17 .
In this paper, we apply this formula to a variety of 1d and 2d models belonging to the classes BDI, DIII and D. We show that we can successfully chart out the topological phase diagrams. Studying these examples also enables us to explicitly observe the correspondence between these complex momentum solutions in the Fourier space, and the localized Majorana fermion wavefunctions in the position space. Finally, we also corroborate the fact that for systems with a chiral symmetry, these solutions can be identified with the so-called "Exceptional Points" (EP 's) [18][19][20][21][22][23][24] , where two or more eigenvalues of the complexified Hamiltonian coalesce. EP 's are singular points at which the norm of at least one eigenvector vanishes, when certain real parameters appearing in the Hamiltonian are continued to complex values, and the complexified Hamiltonian becomes non-diagonalizable. The concept of EP 's is similar to that of a degeneracy point, but with the important difference that all the energy eigenvectors cannot be made orthogonal to each other. In previous works, EP 's have been used [25][26][27][28] to describe topological phases of matter for 1d topological superconductors/superfluids. The paper is organized as follows: In Sec. II, we review the results obtained earlier 15 for counting the number (n) of Majorana zero modes bound to defects, based on the bulk-edge correspondence. In Sec. III, we consider some 1d and 2d models in the class BDI and apply the EP formalism to count n. Sec. IV is devoted to the study of edge states for Hamiltonians in class DIII, where we illustrate the applicability of EP solutions as the chiral symmetry exists. In Sec. V, we discuss some systems in the class D and conclude that EP 's cannot be related to the Majorana fermion wavefunctions for such Hamiltonians, as chiral symmetry is broken. We conclude with a summary and outlook in Sec. VI.

II. COUNTING FORMULA FOR THE MAJORANA ZERO MODES
In this section, we review the connection 15 between the complex momentum solutions of det (H BdG (k)) = 0, and the Majorana fermion wavefunctions in the position space associated with a defect or edge.
We consider a topological defect embedded in (or at the boundary of) a d-dimensional topological superconductor. Let m be the dimensions of the defect, parametrized by the Cartesian coordinates r ⊥ = (r 1 , . . . , r d−m ) and r = (r d−m+1 , . . . , r d ), located at r ⊥ = 0. Let k ⊥ = k ⊥Ω = (k 1 , . . . , k d−m ) and k = (k d−m+1 , . . . , k d ) be the corresponding conjugate momenta, where k ⊥ = |k ⊥ | and Ω is the unit vector when written in spherical coordinates.
For a generic H BdG , let k j A and k j B (j = 1, . . . , Q) be the two sets of complex k ⊥ -solutions for det (H BdG (k)) = 0, related by { (k j A )} = −{ (k j B )}, after k ⊥ has been analytically continued to the complex plane. Assuming the Majorana wavefunction to be of the form ∼ exp (−z |r ⊥ |) in the bulk, the correspondence i k ⊥ ↔ −z has been established 15 . At a topological phase transition point, one or more of the k j A/B 's go through zero. When k j A/B changes sign at a topological phase transition point, the position space wavefunction of the corresponding Majorana fermion changes from exponentially decaying to exponentially diverging or vice versa. If the former happens, the Majorana fermion ceases to exist. A new Majorana zero mode appears in the latter case. The count (n) for the Majorana fermions for a defect is thus captured by the function where ({λ i }, k ,Ω) are the parameters appearing in the are their values at any point in the non-topological phase.
In addition to having the PHS and time-reversal symmetry (TRS), if there is a chiral symmetry operator O which anticommutes with the Hamiltonian, the latter takes the form in the momentum space, for the corresponding bulk system with no defect. On analytically continuing the magnitude k ⊥ ≡ k = |k| to the complex k ⊥ -plane, at least one of the eigenvectors of H chiral (k) collapses to zero norm where det (A(k)) = 0 , or det A † (k) = 0 .
These points are associated with the solutions of EP 's for complex k ⊥ -values where two or more energy levels coalesce. Furthermore, these coalescing eigenvalues have zero magnitude as det (A(k)) = 0 (or det A † (k) = 0) also implies det (H chiral (k)) = 0. H chiral (k) becomes nondiagonalizable, as in the complex k ⊥ -plane, det (A(k)) = 0 det A † (k) = 0 (or vice versa). However, at the physical phase transition points, the imaginary parts of one or more solutions vanish, and det (A(k)) = det A † (k) = 0 for those solutions, making H chiral (k) once again diagonalizable and marking the disappearance of the corresponding EP 's.
Because it satisfies Eq. (3), each EP solution corresponds to a Majorana fermion of a definite chirality with respect to O. If A † (k) = A T (−k) holds, then the two sets of EP 's are related by {k j A } = −{k j B }, one set corresponding to the the solutions obtained from one of the two off-diagonal blocks. In such cases, the pairs of the Majorana fermion wavefunctions are of opposite chiralities.

III. EP FORMALISM FOR THE BDI CLASS
In this section, we consider some 1d and 2d spinless models in the BDI class, which can support multiple Majorana fermions at any end of an open chain. For systems in this class, there exists a chiral symmetry operator O, such that H BdG can be rotated to the form H chiral in Eq. (2).
After reviewing the transfer matrix scheme to find Majorana fermion solutions localized at an edge, we show how EP solutions in the complex k ⊥ -plane can be used to count the number of Majorana zero modes in a given topological phase. We also emphasize on the connection of these EP solutions with the position space wavefunctions calculated in the real space lattice with open ends.

A. Transfer matrix approach
Kitaev 2 suggested the model of a 1d p-wave superconducting chain, which can support Majorana zero modes at the two ends. For a finite and open chain with N sites, the Hamiltonian takes the form where µ is the chemical potential, w and ∆ are the nearest-neighbour hopping amplitude and superconducting gap respectively. The pair of fermionic annihilation and creation operators, c j and c † j , describe the lattice site j, and obey the usual anticommutation relations {c j , c j } = 0 and {c j , c † j } = δ jj . The Majorana mode structure of the wire can be understood better by rewriting the above Hamiltonian in terms of the Majorana operators satisfying Then the Hamiltonian reduces to This chain can support one Majorana bound state (MBS) at an edge for appropriate values of the parameters. More recently, a variation of the model was considered with next-nearest-neighbour hopping and pairing amplitudes 29 . A general version of such longer-ranged interactions with all possible hoppings and pairings was studied 30,31 with the Hamiltonian where the J ±r 's are real parameters, and 0 < q < N . These models can support multiple MBSs at an edge. If we impose periodic boundary conditions (PBC's), the Hamiltonian can be diagonalized by a Bogoliubov transformation: where the anticommuting fermion operators c † k , c k are suitable linear combinations in the momentum space of the original c j , c † j fermion operators. The energy eigenvalues are given by We now review the transfer matrix approach 29,31,32 to identify the number of MBSs at each end of the chain for this model. The transfer matrix can be obtained from the Heisenberg equations of motion for the Majorana operators in Eq. (9): Assuming the time-dependence to be of the form a j = A j e −iE l t and b j = B j e −iE l t , the E l = 0 (zero energy modes) are given by the recursion relation of the amplitudes: Clearly, it will suffice to solve one set of the recursive equations to obtain the solutions for both. Assuming An MBS can exist if we have a normalizable solution, i.e., if |λ A | < 1 or |λ B | < 1, if the solution is to be localized at the left end. Similarly, for a mode to be localized at the right end of the chain, we must have |λ A | > 1 or |λ B | > 1. Depending on the number of constraint equations (or boundary conditions on the amplitudes), one has to determine the number of independent MBSs at each end of the chain.

B. Relation of the EP formalism with the transfer matrix approach
Let us apply the EP formalism 15,28 to the Hamiltonian in Eq. (11). First we rotate it to the off-diagonal form The EP 's where either A l (k) or B l (k) vanishes, are given by the solutions Comparing Eqs. (16), (18) and (19), it is easy to see that the solutions for EP 's in the complex k-plane for the PBC's correspond to the MBS solutions for the open boundary conditions (OBC's). Since a sign change of k Al/Bl indicates a topological phase transition, by which we move from a phase where an MBS can exist to one where that particular zero mode gets destroyed. This is related to the fact that k Al/Bl 's are related to the exponential decay of the MBS position space wavefunctions localized at one end of the open chain.
and all other J r 's to be zero, we can get a system supporting upto four Majorana zero modes at each end of the chain. The phase diagram obtained using Eq. (1) is shown in Fig. 1.
Instead, for the parameters , and all other J r 's set to zero, we get a system having three EP 's for either  We should note another important point: if there are Q EP solutions for either A l (k) = 0 or B l (k) = 0, clearly there are 2 Q solutions in total. But for counting the zero modes in Eq. (1), we should consider only one set, where the two sets obey the relatioñ As we have already seen, these two sets correspond to the wavefunctions of the MBSs at the two opposite ends. Evidently, the MBSs exist in pairs at the two ends and the topological phase is characterized by their number at each individual end.

C. Single-channel ferromagnetic nanowire
The 1d Hamiltonian for a ferromagnetic nanowire embedded on Pb superconductor 33 with a single spatial channel (i.e. no transverse hopping) is given by Here k is the 1d crystal momentum, Ψ k is the fourcomponent Nambu spinor defined in the particle-hole (τ ) and spin (σ) spaces, and V is the Zeeman field which can be induced by ferromagnetism. Also, ∆ s and ∆ p are proximity-induced s-wave and p-wave superconducting pairing potentials respectively, with d determining the relative magnitudes of the components of the p-wave superconducting order parameter ∆ αβ (α, β =↑, ↓). In our calculations, we use d = (1, 0, 0) and V = (0, 0, V ). This Hamiltonian then belongs to the BDI class with the chiral symmetry operator given by O = σ x τ y .
The eigenvalues of the Hamiltonian are given by: A level crossing can occur if either E 1 (k) = 0 orẽ = 0. However, for a finite V and ∆ s , the latter is impossible. Hence, a level crossing takes place when E 1 (k) = 0 for k = 0 or π for the appropriate values of the parameters, which also indicates that this corresponds to the appearance of zero energy modes. A generic complex value of k corresponding to E 1 (k) = 0 can be obtained by solving We can rotate the Hamiltonian in Eq. (22) to the chiral basis, where it takes the form If either det(h N 1 u (k)) = 0 or det(h N 1 l (k)) = 0 for a complex k-value, this leads to the vanishing of the norm of one of the four eigenvectors of H chi N 1 (k), signalling the existence of an EP for that value of k. At an EP , H chi N 1 (k) is thus non-diagonalizable.
The solutions for the EP 's are given by either where (s 1 = ±1, s 2 = ±1). Clearly, k = k u/l s1,s2 also solves Eq. (24), which corresponds to two coinciding zero energy solutions (where two levels coalese 22 for a complex kvalue).
The plots of the energy bands, k u s1,s2 , and f (µ) have been shown in Fig. 3 , using the values ∆ s = ∆ p = 0.1 t andṼ = 1.5 t. Now let us try to understand the existence of the EP 's throughout a given a topological phase and their disappearance right at the phase transition points, the latter being tied to the sign change of the k u/l s1,s2 's. The Hamiltonian in Eq. (25), when written in position space, gives the following equations for the Majorana zero modes, ψ + = (u + , 0) T and ψ − = (0, u − ) T (with chirality +1 and −1 respectively): Here we have assumed a continuum for an open wire and set t = 1. For ψ − , let us assume the trial solution The complex z r 's must satisfy the quartic equation whereas for small k, from Eq. (27), we get indicating correspondence i k ↔ −z r . The magnitude of z r will determine the admissible MBS solutions subject to OBC's (as analyzed in an earlier work 34 ), just as in the transfer matrix analysis for the 1d spinless lattice case. Hence, here also we have been able to establish the relation between the existence of EP 's in the complex k-plane (for the periodic Hamiltonian) and the localized Majorana zero modes at the ends of an open chain.

D. Two-channel time-reversal-symmetric nanowire system
MBSs in a two-channel TRS nanowire proximitycoupled to an s-wave superconductor have been recently studied 35 . The low-energy model for the lowest bands of the system is described by the effective 1d 4 × 4 BdG Hamiltonian: where p c is the momentum when the gap closes, B is a magnetic field for the Zeeman term, and (v,μ) are effective parameters. We have set p c = 2vm for our calculations. Since this nanowire system belongs to the BDI class when B is perpendicular to the spin-orbit-coupling direction, we will take B = (B, 0, 0) in our analysis. Then the chiral symmetry operator is given by O = σ z τ y . The eigenvalues of the Hamiltonian are given by: We can have two levels coalescing if E 1 (k) = 0 for a complex k-value obtained by solving As before, we rotate the Hamiltonian in Eq. (32) to the chiral basis, where it takes the form The solutions for the EP 's are then given by either where (s 1 = ±1, s 2 = ±1). Clearly, k = k u/l s1,s2 also solves Eq. (34) and hence correspond to the coalescing of two energy levels at the zero value in the complex k-plane.
Choosing v = 1 and m = 1 2 v 2 , the energy bands for B = 0 and B = 3, and the contourplot for f (μ, B) (defined in Eq. (1) have been shown in Fig. 4. Once again we find that f (μ, B) gives the correct topological phase diagram in Fig. 4(c). Needless to add that here also exp i k  In this subsection, we consider the EP -formalism for a 2d-lattice Hamiltonian in the class BDI. The Kitaev honeycomb model 36 can be mapped onto free spinless fermions with p-wave pairing on a honeycomb lattice, using the Jordan-Wigner transformation. The solutions for the edge modes for a semi-infinite lattice have been studied earlier [37][38][39][40] . The momentum space Hamiltonian in terms of the Majorana operators is with the eigenvalues We will consider two kinds of edges 37,38 , namely, zigzag and armchair, which can support Majorana fermions. We will find the phase diagram using the EP 's corresponding to these edges setting either A(k) = 0 or B(k) = 0, after complexifying the momentum component perpendicular to the edge. For this 2d case, f in Eq. (1) is a function of (J 1 , J 2 , J 3 , k ), where k is the momentum along the 1d edge being considered.
One can have a zigzag edge in the y-direction, according to the convention of Nakada et al 37 , so that we will complexify k ⊥ = k x , and k = k y will be one of the parameters determining the topological phase transition points. The solution for B(k x = k ⊥ , k y = k ) = 0 is given by For the armchair edge 37 in the x-direction on the top of the lattice, we will complexify k ⊥ = k y , and k = k x will be now one of the parameters determining the topological phase transition points. The two EP 's for A(k x = k , k y = k ⊥ ) = 0 are given by Eqs. (40) and (41) This system may be realized in a Rashba wire that is proximity-coupled to a nodeless s ± wave superconductor. The energy eigenvalues are given by: whose plots are shown in Fig. 6(a) for ∆ = 0.1 t and λ R = 2 t, as µ/t is varied along the horizontal axis.
Observing that a chiral symmetry operator O = σ 0 τ y exists in the presence of M z , we rotate the Hamiltonian in Eq. (42) to the chiral basis, where it takes the form The solutions for the EP 's are then given by either or det(h l m1 (k)) = 0 ⇒ i ∆ cos(k) + ξ m1 (k) = s 1 λ R sin(k) where (s 1 = ±1, s 2 = ±1). For this model, we note that the two different sets of EP solutions, related by { (k)}| set=A = −{ (k)}| set=B , are obtained from Eqs. (46) and (47) when we set s 1 = 1 and s 1 = −1 respectively. This is related to the fact that h u m1 (k) However, we have argued before that for Eq. (1) to work, we must take all the EP solutions from one of the sets related by a negative sign of (k). Using either (k m1u +1,s2 , k m1l +1,s2 ) or (k m1u −1,s2 , k m1l −1,s2 ) (but not both), Fig. 6(b) gives the correct topological phase diagram in the µ/t − λ R /t plane, for ∆ = 0.1 t. We clearly see that there exists phases with a pair of MBSs, which correspond to one Kramers doublet (MKP).

B. 2d model
The 1d model of a Rashba semiconductor combined with a nodeless s ± wave superconductor can be easily generalized to a 2d system, described by the Hamiltonian 42 where ∆ m (k) is the s ± wave singlet pairing potential that switches its sign between the centre (0, 0) and the corner (π, π) of the 2d Brillouin zone, when 0 < |∆ 0 | < 4∆ 1 .

V. BROKEN TIME REVERSAL SYMMETRY: CLASS D
In this section, we consider 1d and 2d Hamiltonians in the symmetry class D, where the TRS is broken. We will see that the EP formalism in the complex k ⊥ -plane is not applicable for such systems.
A. 1d spinless model We examine the spinless model described by the Hamiltonian which looks similar to H K in Eq. (4), but with the important difference that the TRS is broken by the fact that w and ∆ can be complex numbers 31 . Without any loss of generality, we can choose ∆ to be real and encode the entire phase-difference (φ) between ∆ and w 0 by writing w = w 0 e iφ , where w 0 is real and positive. With PBC's, one can write the corresponding BdG Hamiltonian in the momentum space as: The energy eigenvalues are given by: e t = (2 w 0 cos(φ) cos(k) + µ) 2 + 4 ∆ 2 sin 2 (k) .
Hence it follows that two levels become degenerate wheñ e t = 0 ⇒ 2 w 0 cos(φ) cos(k) + µ = ± 2 i ∆ sin(k) . (61) When extended to the complex k-space, there can be EP 's where the norm of an eigenvector vanishes. For convenience, we rotate h D1 (k) to find the points (EP 's) where the Hamiltonian becomes non-diagonalizable: The EP 's are given by A D1 (k) = 0 or B D1 (k) = 0. At such points, two levels coalesce for a complex value of k satisfying Eq. (61).
However, we immediately observe that these EP 's do not correspond to zero energy modes, which appear when (2 w 0 cos(φ) cos(k) + µ) 2 + 4 ∆ 2 sin 2 (k) is satisfied. Although the EP description no longer applies now to the existence of MBSs, we can find the complex k-values satisfying det (h D1,od (k)) = 0 We can solve for the k-values either with the "+" or the "−" sign on the RHS (but not both), and plug in the roots of that equation into the formula in Eq. (1). The function f (µ, ∆) will still give the number of MBSs in a given topological phase. Once again we emphasize that for counting the zero modes in Eq.
(1), we should include only one of the two sets of roots related by a sign change of (k), as these two sets correspond to the wavefunctions of the pair of MBSs at the two opposite edges. We have shown the plots of f (µ, ∆) in Fig. 8 for four different values of φ. It indeed captures the correct Z 2 topological invariant (n = 0, or 1). No zero mode exists, i.e., the system is entirely gapped in certain regions in the φ − µ plane where (k) vanishes, as (k) is related to the exponential part of the MBS wavefunction in the real space.
B. 2d spinless model The following Hamiltonian gives a model of a p + ip wave superconductor on a square lattice 43 : h D2 (k) = (2 t x cos(k x ) + 2 t y cos(k y ) − µ) τ z + d x sin(k x ) τ x + d y sin(k y ) τ y , where (t x , t y ) are the hopping strengths and (d x , d y ) are the pairing amplitudes along the (x, y)-directions, and µ is the chemical potential. The energy eigenvalues, given by E(k) = ± (2 t x cos(k x ) + 2 t y cos(k y ) − µ) 2 + d 2 x sin 2 (k x ) + d 2 y sin 2 (k y ) are plotted in Fig. 9(a) for t x = t y = d x = d y = 1, as µ is varied along the horizontal axis. We consider edges parallel to the y-axis, so that k = k y and k ⊥ = k x . Complexifying k x , we can compute the number of non-chiral Majorana fermions propagating along these edges with momenta k y = 0 and k y = π, using Eq. (1). The complex k x -values satisfying det (h D2 (k)) = 0, for k y = 0 and k y = π, are given by respectively. Here (s 1 = ±1, s 2 = ±1), and we need to use either s 1 = 1 or s 1 = −1 (but not both) in Eq. (1), to obtain the phase diagrams shown in Figs. 9(b) and 9(c).

VI. CONCLUSION
We have established the relation of the EP solutions for complexified momenta to the Majorana fermion wavefunctions bound to a topological defect in a system with a chiral symmetry, by studying some explicit examples in 1d and 2d. These models include both spinless and spinful cases. We have shown that such EP solutions cannot exist for systems in class D, where the chiral symmetry is absent. The generic formula, which was proposed earlier 15 to count the number of Majorana zero modes in arbitrary dimensions, has been demonstrated to chart out the desired topological phase diagrams for the wide variety of systems we have considered. The detailed study of these models also helps us illustrate how one distinct set of complex k ⊥ -solutions for det (H BdG (k)) = 0, related by { (k ⊥ )}| set=A = −{ (k ⊥ )}| set=B , has to be used while using our formula. For a system with or without a chiral symmetry, the imaginary parts of these solutions in the complexified k ⊥ -plane are related to the exponential decay of the Majorana fermion wavefunctions in the bulk in the position space. Hence, the imaginary parts of n of the solutions in one set undergoes a change of sign across a topological phase transition point, if the number of Majorana zero modes at a defect changes byñ.