3. Derivation of the counting formula

DOI: 10.5488/CMP.19.33703 http://www.icmp.lviv.ua/journal

Counting Majorana bound states using complex momenta

I. Mandal

Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo ON N2L 2Y5, Canada

Received February 17, 2016, in final form April 14, 2016

Recently, the connection between Majorana fermions bound to the defects in arbitrary dimensions, and com- plex momentum roots of the vanishing determinant of the corresponding bulk Bogoliubov–de Gennes (BdG) Hamiltonian, has been established (EPL, 2015,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 correspon- dence 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 symme- try, these solutions are the so-called “exceptional points”, where two or more eigenvalues of the complexified Hamiltonian coalesce.

Key words:exceptional points, Majorana fermions, BDI, DIII, D, counting PACS:73.20.-r, 74.78.Na, 03.65.Vf

1. Introduction

Topological superconductors [1] are systems which can provide the condensed matter version of Ma- jorana fermions, because 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 localized zero-energy 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 di- mensional spinlessp-wave superconducting systems [2, 5], as well as other systems involving various heterostructures with proximity-induced superconductivity which are topologically similar to them [6–

12].

Non-interacting Hamiltonians for gapped topological insulators and topological superconductors, in arbitrary spatial dimensions, can be classified into ten topological symmetry classes [13–15], character- ized by certain topological invariants. Moreover, there exists a unified framework for classifying topo- logical defects in insulators and superconductors [16], which follows from the bulk-boundary correspon- dence and identification of the protected gapless fermion excitations with topological invariants char- acterizing the defect. 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 [13] into three classes:

BDI, DIII and D.

In our earlier work [17], we have explored the connection between the complex momentum solu- tions of the determinant of a bulk BdG Hamiltonian (HBdG) 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 familiarZ^andZ2topological invariants [13, 14, 18] and other counting schemes [19–22].

In this paper, we prove this formula and apply it 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. Study- ing these examples also enables us to explicitly observe the correspondence between these complex mo- mentum 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) [23–29], 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 [30–34] to describe topological phases of matter for 1d topological superconductors/superfluids.

The paper is organized as follows: in section 2, we review the results obtained earlier [17] for counting the number (n) of Majorana zero modes bound to defects, based on the bulk-edge correspondence. In section 3, we provide a proof of the counting formula. In section 4, we consider some 1d and 2d models in the class BDI and apply the EP formalism to countn. Section 5 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 section 6, we discuss some systems in the class D and conclude that EP’s cannot be related to the Majorana fermion wavefunctions for such Hamiltonians, because chiral symmetry is broken. We conclude with a summary and outlook in section 7. In appendix A, we provide a simple example to show how one should choose the correct EP solutions such that their imaginary parts are continuous functions in the parameter space in order to evaluate our counting formula.

2. Counting formula for the Majorana zero modes

In this section, we review the connection [17] between the complex momentum solutions of det[HBdG(_k)]=0, and the Majorana fermion wavefunctions in the position space associated with a de- fect or edge.

We consider a topological defect embedded in (or at the boundary of ) a d-dimensional topologi- cal superconductor. Letmbe the dimensions of the defect, parametrized by the Cartesian coordinates r_⊥=(r1, . . . ,r_d−m)^andr_∥=(r_d−m+1, . . . ,r_d), located atr_⊥=0^{. Let}k_⊥=k_⊥Ωˆ =(k1, . . . ,k_d−m)^andk_∥= (k_d−m+1, . . . ,k_d)be the corresponding conjugate momenta, wherek_⊥= |k_⊥|^and Ωˆ is the unit vector when written in spherical coordinates.

For a generic HBdG, let k_A^{j and} k_B^{j (}j = 1, . . . ,Q) be the two sets of complex k_⊥-solutions for det[HBdG(_k)]=0, related by{Im(k_A^j)}= −{Im(k_B^j)}^{, after}k_⊥has been analytically continued to the complex plane. One should be careful to choose solutions such that their imaginary parts are continuous functions of the parameter(s) which tune(s) through the transition, and the solutions in one set are related to the other by changing the sign of their imaginary parts throughout. This point has been illustrated by an example in appendix A. Assuming the Majorana wavefunction to be of the form∼exp (−z|r⊥|)^{in the} bulk, the correspondenceik_⊥↔ −zhas been established [17]. At a topological phase transition point, one or more of theIm(k_A/B^j )’s go through zero. WhenIm(k_A/B^j )changes sign at a topological phase transition point, the position space wavefunction of the corresponding Majorana fermion changes from exponen- tially 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 captured by the function

f({λi},_k∥, ˆΩ)=1 2

¯

¯

¯

¯

Q

X

j=1

µ sign

n Imh

k_A/B^j ¡

{λi},_k∥, ˆΩ¢io

−^signn Imh

k_A/B^j ¡

{λ⁰_i},_k⁰_∥, ˆΩ⁰¢io¶ ¯

¯

¯

¯

, ^(2.1)

where({λi},_k∥, ˆΩ)are the parameters appearing in the expressions fork_A/B^{j , and}({λ⁰_i},_k⁰_∥, ˆΩ⁰)^{are their} values at any point in the non-topological phase.

If there is a chiral symmetry operatorO which anticommutes with the Hamiltonian, the latter takes the form

Hchiral(_k)=

µ 0 A(k) A^†(_k) 0

¶

, ^(2.2)

in the momentum space, for the corresponding bulk system with no defect. On analytically continuing the magnitudek_⊥≡k= |k|to the complexk_⊥-plane, at least one of the eigenvectors ofHchiral(_k)^collapses to zero norm where

det [A(_k)]=0 ^or det£ A^†(_k)¤

=0. ^(2.3)

These points are associated with the solutions of EP’s for complexk_⊥-values where two or more energy levels coalesce. Furthermore, these coalescing eigenvalues have zero magnitude sincedet[A(_k)]=0^(or det[A^†(_k)]=0) also impliesdet[Hchiral(_k)]=0^.Hchiral(_k)becomes non-diagonalizable, as in the com- plexk_⊥^-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, anddet[A(_k)]=det[A^†(_k)]=0for those so- lutions, makingHchiral(_k)once again diagonalizable and marking the disappearance of the corresponding EP’s.

Since it satisfies equation (2.3), each EP solution corresponds to a Majorana fermion of a definite chirality with respect toO^{. If}A^†(_k)=A^T(−k)holds, then the two sets of EP’s are related by{k^j_A}= −{k_B^j}^, one set corresponding to 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.

3. Derivation of the counting formula

A simple derivation of the counting formula in equation (2.1) can be motivated as follows:

1. Let us consider one of the solutions given by j=1. In the non-topological phase, say phase “0^”, k_A¹({λ⁰_i},_k⁰_∥, ˆΩ⁰)gives no Majorana zero mode and hence does not give rise to any decaying mode localized at a defect. On the other hand, in a topological phase, say phase “t”, with a Majorana wavefunction ∼exp[−|Im(k¹_A)|r_⊥]^,k¹_A({λi},_k∥, ˆΩ)¯

¯_{phase t} localized atr_⊥=0and zero atr_⊥= ∞^, should now give rise to an admissible decaying zero mode solution. This implies that there is a change in sign ofIm(k_A¹)^from−1^to+1when we jump from phase “0” to phase “t”.

2. Majorana zero modes must occur in pairs, though they might be localized far apart. Hence, if k_A¹¯

¯_{phase C}corresponds to a Majorana mode localized atr_⊥=0^{, then}k_B¹¯

¯_{phase t}must correspond to one localized atr_⊥= ∞^{, where}k¹_B=(k¹_A)^∗. Hence, whether or not we are in the topological phase “t”

is captured by the functionf₁=¹₂¯

¯^sign

©Im[k¹_A/B({λi},_k∥, ˆΩ)]ª

−^sign©

Im[k_A/B¹ ({λ⁰_i},_k⁰

∥, ˆΩ⁰)]ª¯

¯^taking the value1^{or zero.}

3. From the above discussion, it may seem that the counting formula should be given by 1

2

Q

X

j=1

¯

¯

¯^sign n

Imh k_A/B^j ¡

{λi},_k∥, ˆΩ¢io

−^signn Imh

k_A/B^j ¡

{λ⁰_i},_k⁰_∥, ˆΩ⁰¢io ¯

¯

¯.

However, this is not quite correct. To understand this, let us consider the scenario when at least two of the solutions, sayk_A^1andk_A²are such that sign[Im(k_A²)]¯

¯_phase0= −^sign[Im(k¹_A)]¯

¯_phase0. This implies that in the trivial phase, the wavefunction given byc₁exp(ik¹_A¯

¯_phase0r_⊥)+c₂exp(ik_A²¯

¯_phase0r_⊥)^is inadmissible for not being capable of satisfying the boundary conditions — the only solution is c₁=c₂=0. In another topological phase, say “˜t ”, let sign[Im(k_A¹)]¯

¯_phase˜t= −^sign[Im(k¹_A)]¯

¯_phase0^and sign[Im(k_A²)]¯

¯_phase˜t= −^sign[Im(k²_A)]¯

¯_phase0. This means that bothIm(k_A¹)^andIm(k²_A)change sign when we jump from phase “0” to phase “˜t ”. However, they still should not give any Majorana zero mode in the phase “˜t ”, becausec˜₁exp(ik_A¹¯

¯_phase˜tr_⊥)+c˜₂exp(ik_A²¯

¯_phase˜tr_⊥)cannot satisfy the boundary conditions. So, the correct formula is given by equation (2.1).

4. 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 operatorO, such thatHBdGcan be rotated to the formHchiralin equation (2.2).

After reviewing the transfer matrix scheme to find Majorana fermion solutions localized at an edge, we show how EP solutions in the complexk_⊥-plane can be used to count the number of Majorana zero modes in a given topological phase. We also make emphasis on the connection of these EP solutions with the position space wavefunctions calculated in the real space lattice with open ends.

4.1. Transfer matrix approach

Kitaev [2] suggested the model of a 1dp-wave superconducting chain, which can support Majorana zero modes at the two ends. For a finite and open chain withN sites, the Hamiltonian takes the form

HK= −

N

X

j=1

µ µ

c^†_jcj−1 2

¶ +

N−1X

j=1

³

−w c^†_jcj+1+∆cjcj+1+^h.c.´

, ^(4.1)

whereµis the chemical potential,w^and∆are the nearest-neighbour hopping amplitude and supercon- ducting gap, respectively. The pair of fermionic annihilation and creation operators,cj andc^†_j^{, describe} the lattice sitej, and obey the usual anticommutation relations{c_j,c⁰_j}=0^and{c_j,c^†

j⁰}=δj j⁰. The Majo- rana mode structure of the wire can be better understood by rewriting the above Hamiltonian in terms of the Majorana operators

a_j=c^†_j+c_j, b_j= −i³ c^†_j−c_j´

, ^(4.2)

satisfying

a_j=a^†_j, b_j=b^†_j, {a_j,b_j0}=0 , {a_j,a_j0}={b_j,b_j0}=2δj j⁰. Then, the Hamiltonian reduces to

HK= −i 2

N

X

j=1

µajbj−i 2

N−1X

j=1

£(w−∆)ajbj+1−(w+∆)bjaj+1¤

. ^(4.3)

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 hop- ping and pairing amplitudes [35]. A general version of such longer-ranged interactions with all possible hoppings and pairings was studied [36, 37] with the Hamiltonian

H_l= −i 2

N

X

j=1

µa_jb_j−i

q

X

r=1 N−q

X

j=1

£J_−ra_jb_j+r+J_ra_j+rb_j¤

, ^(4.4)

where theJ_±r’s are real parameters, and0<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:

H_l= −X

k

³

c_k^† c_−k ´ h_l(k)

µ c_k c_−k^†

¶ ,

h_l(k)= −2

q

X

r=−q

µ J_rcos (kr) −iJ_r sin (kr) iJrsin (kr) −Jr cos (kr)

¶

, J₀= −µ

2, ^(4.5)

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

E_l(k)= ±2 s

hX

r

J_r cos (kr) i2

+ hX

r

J_r sin (kr) i2

. ^(4.6)

We now review the transfer matrix approach [35–38] 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 equation (4.4):

2 ida_j dt = −i

q

X

r=−q

J₋rbj+r, 2 idb_j dt =i

q

X

r=−q

Jraj+r. ^(4.7)

Assuming the time-dependence to be of the forma_j=A_je^−iEltandb_j=B_je^{−iElt, the}E_l=0(zero energy modes) are given by the recursion relation of the amplitudes:

q

X

r=−q

J_−rb_j+r=0 ,

q

X

r=−q

J_ra_j+r=0 . ^(4.8)

Clearly, it will suffice to solve one set of the recursive equations to obtain the solutions for both. Assuming A_j=λA^j andB_j=λB^j, we get the polynomial equations

q

X

r=−q

J_rλ^q+rA =0 ,

q

X

r=−q

J_−rλB^q+r=0 . ^(4.9)

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-hand 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 should determine the number of independent MBSs at each end of the chain.

4.2. Relation of the EP formalism with the transfer matrix approach

Let us apply the EP formalism [17, 34] to the Hamiltonian in equation (4.5). First we rotate it to the off-diagonal form

h_l,od(k)=U^†

l h_l(k)U_l=

µ 0 A_l(k) B_l(k) 0

¶

, U_l= i

p2

µ −1 −1

−1 1

¶ ,

A_l(k)= −2

q

X

r=−q

£J_rcos (kr)+iJ_rsin (kr)¤

, B_l(k)= −2

q

X

r=−q

£J_rcos (kr)−iJ_r sin (kr)¤

. ^(4.10)

The EP’s where eitherA_l(k)^orB_l(k)vanishes, are given by the solutions q

X

r=−q

Jrλ˜A^q+l^r=0, ^where λ˜^Al=exp (ikAl) , ^(4.11)

q

X

r=−q

J_−rλ˜B^q+l^r=0, ^where λ˜^Bl=exp (ikBl) . ^(4.12)

Comparing equations (4.9), (4.11) and (4.12), 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

|λ˜^Al/^Bl| <1 ⇒ Im (kAl/^Bl)>0 ⇔ |λ^{A / B}| <1 , ^(4.13)

|λ˜^Al/^Bl| >1 ⇒ Im (kAl/^Bl)<0 ⇔ |λ^{A / B}| >1 , ^(4.14) a sign change ofIm (kAl/^Bl)indicates a topological phase transition, by which we move from a phase where an MBS can exist to the one where that particular zero mode gets destroyed. This is related to the fact thatIm (kAl/^Bl)’s are related to the exponential decay of the MBS position space wavefunctions localized at one end of the open chain.

Choosing J₀= −^µ₂,J₁=J₂=^1+∆₂ ,J₋₁=J₋₂= ^1−∆₂ and all other J_r’s to be zero, we can get a system supporting up to four Majorana zero modes at each end of the chain. The phase diagram obtained using equation (2.1) is shown in figure 1.

Figure 1. (Color online) The topological phase diagram of the Hamiltonian described by equation (4.4), withJ0= −^µ₂,J1=J2=^1+∆₂ ,J₋1=J₋2=^1−∆₂ , and all otherJr’s set to zero. Here,nlabels the number of Majorana zero modes at each end of the chain, as captured by the functionf(µ,∆)defined in equa- tion (2.1).

Instead, for the parameters J0= −^M₂cos(φ2),J1= −^J₂cos(φ1),J₋1= −₂^Jsin(φ1), and all otherJr’s set to zero, we get a system having three EP’s for eitherA_l(k)=0^orB_l(k)=0. For this model, up to two Ma- jorana zero modes can appear at an edge. The phase diagrams forJ/M=0.625^andJ/M=1.3^{, obtained} using equation (2.1), are shown in figure 2.

We should note another important point: if there areQEP solutions for eitherA_l(k)=0^orB_l(k)=0^, clearly there are2Qsolutions in total. However, for counting the zero modes in equation (2.1), we should consider only one set, where the two sets obey the relation

λ˜^Al=1/ ˜λ^Bl or kAl= −kBl. ^(4.15)

(a) (b)

Figure 2. (Color online) Panels (a) and (b) show the topological phase diagram of the Hamiltonian de- scribed by equation (4.4), withJ₀= −^M₂cos(φ2),J₁= −₂^Jcos(φ1),J₋₁= −₂^Jsin(φ1), and all otherJr’s set to zero. Here,nlabels the number of Majorana zero modes at each end of the chain, as captured by the functionf(µ,∆)defined in equation (2.1).

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.

4.3. Single-channel ferromagnetic nanowire

The 1d Hamiltonian for a ferromagnetic nanowire embedded on Pb superconductor [39] with a single spatial channel (i.e., no transverse hopping) is given by

H_N1=X

k

Ψ^†_kh_N1(k)Ψk, Ψk=(c_k↑,c_k↓,c^†

−k↓,−c^†

−k↑)^T, h_N1(k)=ξ(k)σ0τz+£

∆sσ0+∆psin(k)d·σ¤

τx+V·στ0, ξ(k)= −2tcos(k)−µ. ^(4.16) Here,k is the 1d crystal momentum,Ψk is the four-component Nambu spinor defined in the particle- hole(τ)^{and spin}(σ)spaces, andVis the Zeeman field which can be induced by ferromagnetism. Also,

∆s and∆p are proximity-induceds^{-wave and}p-wave superconducting pairing potentials, respectively, withddetermining the relative magnitudes of the components of the p-wave superconducting order parameter∆αβ (α,β=↑,↓). In our calculations, we used=(1, 0, 0)^and V=(0, 0,V). This Hamiltonian belongs to the BDI class with the chiral symmetry operator given byO=σxτy.

The eigenvalues of the Hamiltonian are given by:

E₁(k)= ±

qξ²(k)+V²+∆²s+∆²psin²(k)−e˜₁, E₂(k)= ±

qξ²(k)+V²+∆²s+∆²psin²(k)+e˜₁,

e˜1=2 q

V²£

∆²s+ξ²(k)¤

+∆²s∆²psin²(k) . ^(4.17)

A level crossing can occur if eitherE₁(k)=0^orE₂(k)=0. However, for a finiteV ^and∆s, the latter is impossible. Hence, a level crossing takes place whenE₁(k)=0^fork=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 ofkcorresponding toE₁(k)=0can be obtained by solving

hV˜²−ξ²(k)+∆²psin²(k)i2

+4ξ²(k)∆²psin²(k)=0 , V˜= q

V²−∆²s. ^(4.18) We can rotate the Hamiltonian in equation (4.16) to the chiral basis, where it takes the form

H_N1^chi(k)=U₂^†h_N1(k)U₂=

µ 0 h^N1_u (k) h^N1_l (k) 0

¶

, U₂=1 2









−1−i 0 1+i 0

0 1+i 0 −1−i

0 1−i 0 1−i

1−i 0 1−i 0







 ,

h^N1_u (k)=

µ −ξ(k)+i∆psin(k)−V ∆s

∆s ξ(k)−i∆psin(k)−V

¶ ,

h_l^N1(k)=

µ −ξ(k)−i∆psin(k)−V ∆s

∆s ξ(k)+i∆psin(k)−V

¶

. ^(4.19)

If eitherdet[h^N_u¹(k)]=0^ordet[h_l^N1(k)]=0for a complexk-value, this leads to the vanishing of the norm of one of the four eigenvectors ofH_N1^chi(k), signalling the existence of an EP for that value ofk. At an EP, H_N1^chi(k)is thus non-diagonalizable.

The solutions for the EP’s are given by either det£

h_u^N1(k)¤

=0 ⇒ V˜²−ξ²(k)+∆²psin²(k)= −2 iξ(k)∆psin(k)

⇒ k=k^u_s1,s2= −i ln

½·

s₁V˜−µ+s₂ q

¡s₁V˜−µ¢2

−4t²

¸

¡2t+∆p¢₋₁

¾

, ^(4.20)

or

det£ h^N_l ¹(k)¤

=0 ⇒ V˜²−ξ²(k)+∆²psin²(k)=2 iξ(k)∆psin(k)

⇒ k=k_s^l_1,s2= −i ln

½·

s1V˜−µ+s2

q

¡s1V˜−µ¢2

−4t²

¸

¡2t−∆p¢₋1¾

, ^(4.21)

-4 -2 2 4

Μ  t

-4 -2 2

E

H k L

(a)

(b)

-6 -4 -2 2 4 6

μ/t 1

2 f(μ)

(c)

Figure 3. (Color online) Parameters: ∆s=∆p =0.1t^,V˜ =1.5t corresponding to the Hamiltonian in equation (4.16). (a) Energy bandsE1,2(k), given in equation (4.17), have been plotted in blue and red, respectively, as functions ofµ/t. (b) Plots ofIm (k^u_s1,s2)^versusµ/t^{. (c)}f(µ)giving the count of the chiral Majorana zero modes as a function ofµ/t^.

where(s₁= ±1^,s₂= ±1)^{. Clearly,}k=k^u/l_s

1,s2also solves equation (4.18), which corresponds to two coincid- ing zero energy solutions (where two levels coalesce [27] for a complexk^-value).

The plots of the energy bands,Im (k^u_s1,s2)^{, and} f(µ)have been shown in figure 3 , using the values

∆s=∆p=0.1t^andV˜=1.5t^.

Now, let us try to understand the existence of the EP’s throughout a given topological phase and their disappearance right at the phase transition points, the latter being tied to the sign change of the Im (k^u/l_s

1,s₂)’s. The Hamiltonian in equation (4.19), when written in position space, gives the following equations for the Majorana zero modes,ψ+=(u₊, 0)^Tandψ−=(0,u₋)^T(with chirality+1^and−1^{, re-} spectively):

µ ∂²x+µ+∆p∂x−V ∆s

∆s −∂²x−µ−∆p∂x−V

¶ u₊=0 ,

µ ∂²x+µ−∆p∂x−V ∆s

∆s −∂²x−µ+∆p∂x−V

¶

u₋=0 . ^(4.22)

Here, we have assumed a continuum for an open wire and set t=1^{. For}ψ−, let us assume the trial solutionu₋=X

r

exp (−z_rx) Ã u^↑_r

u^↓_r

!

. The complexz_r’s must satisfy the quartic equation

det

µ z²_r+µ−∆pz_r−V ∆s

∆s −z_r²−µ+∆pzr−V

¶

=0 ⇒ ¡

z_r²+µ−∆pz_r¢2

=V˜², ^(4.23)

whereas for smallk, from equation (4.21), we get

¡−k²+µ+i∆pk¢2

=V˜², ^(4.24)

indicating correspondenceik↔ −zr. The magnitude ofzr will determine the admissible MBS solutions subject to OBC’s (as analyzed in an earlier work [40]), just as in the transfer matrix analysis for the 1d spin- less lattice case. Hence, here also we have been able to establish the relation between the existence of EP’s in the complexk-plane (for the periodic Hamiltonian) and the localized Majorana zero modes at the ends of an open chain.

4.4. Two-channel time-reversal-symmetric nanowire system

MBSs in a two-channel TRS nanowire proximity-coupled to ans-wave superconductor have been re- cently studied [41]. The low-energy model for the lowest bands of the system is described by the effective 1d4×4BdG Hamiltonian:

HN2=X

k

Ψ^†_khN2(k)Ψk, Ψk=(c_k↑,c_k↓,c_−k↓^† ,−c_−k↑^† )^T,

hN2=ξ˜(k)σ0τz+v¡

kσz−pcσ0¢

τx+B·στ0, ξ˜(k)= k²

2m−µ˜, ^(4.25)

wherep_cis the momentum when the gap closes,Bis a magnetic field for the Zeeman term, and(v, ˜µ)^are effective parameters. We have setpc=2vmfor our calculations. Since this nanowire system belongs to the BDI class whenBis perpendicular to the spin-orbit-coupling direction, we will takeB=(B, 0, 0)^{in our} analysis. Then, the chiral symmetry operator is given byO=σzτy.

The eigenvalues of the Hamiltonian are given by:

E₁(k)= ± q

B²+k²v²+4m²v⁴+ξ˜²(k)−e˜₂, E₂(k)= ± q

B²+k²v²+4m²v⁴+ξ˜²(k)+e˜₂,

˜ e2=2

q

4m²v⁴¡

B²+k²v²¢

+B²ξ˜²(k) . ^(4.26)

We can have two levels coalescing ifE₁(k)=0for a complexk-value obtained by solving

£B²+k²v²−ξ˜²(k)−4m²v²¤2

+4 ˜ξ²(k)k²v²=0 . ^(4.27) As before, we rotate the Hamiltonian in equation (4.25) to the chiral basis, where it takes the form

H_N2^chi(k)=U₃^†h_N2(k)U₃=

µ 0 h^N2_u (k) h^N2_l (k) 0

¶

, U₃=1 2









0 −1−i 0 1+i

0 1−i 0 1−i

1+i 0 −1−i 0

1−i 0 1−i 0







 ,

h^N2_u (k)=

µ −ξ˜(k)+i¡

kv+2mv²¢

B B −ξ˜(k)+i¡

kv−2mv²¢

¶ ,

h_l^N2(k)=

µ −ξ˜(k)−i¡

kv+2mv²¢

B B −ξ˜(k)−i¡

kv−2mv²¢

¶

. ^(4.28)

The solutions for the EP’s are then given by either det£

h^N2_u (k)¤

=0 ⇒ £ξ˜(k)−ikv¤2

=B²−4m²v⁴

⇒ k=k^u_s1,s2= −i ln µ

s₁ q

m²v²−2mµ˜+2is₂mp

4m²v⁴−B²+imv

¶

, ^(4.29)

or

det£ h^N2_l (k)¤

=0 ⇒ £ξ˜(k)+ikv¤2

=B²−4m²v⁴

⇒ k=k^l_s1,s2= −i ln µ

s1

q

m²v²−2mµ+˜ 2is2mp

4m²v⁴−B²−imv

¶

, ^(4.30)

-10 -5 5 10̎

-5 5

E_1,2HkL

(a)

-6 -4 -2 2 4 6 ̎

-5 5

E_1,2HkL

(b)

n=0 n=2

n=1

n=1

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3

μ˜ /(vpc) B/(vpc)

(c)

Figure 4. (Color online) Parameters:v=1^,m=1/(2v²)^,p_c=2vmcorresponding to the Hamiltonian in equation (4.25). Panels (a) and (b) show the energy bandsE_1,2(k), given in equation (4.26), as functions ofµ˜^{, for}B=0^andB=3, respectively.E_1,2(k)have been plotted in blue and red, respectively. Panel (c) shows the contourplot off(µ)giving the count “n” of the MBSs in theµ/(v p˜ c)−B/(v pc)^plane.

where(s₁= ±1,s₂= ±1)^{. Clearly,}k=k^u/l_s1,s2 also solves equation (4.27) and hence corresponds to the coa- lescing of two energy levels at the zero value in the complexk-plane. Choosingv=1^andm=1/(2v²)^, the energy bands forB=0^andB=3, and the contourplot forf( ˜µ,B)[defined in equation (2.1)] have been shown in figure 4. Once again we find that f( ˜µ,B)gives the correct topological phase diagram in figure 4 (c). Needless to add that here alsoexp(ik^u/l_s

1,s2)’s determine the admissible solutions for the MBS wavefunctions in the position space, at the ends of an open chain.

4.5. Majorana edge modes for the Kitaev honeycomb model

In this subsection, we consider the EP-formalism for a 2d lattice Hamiltonian in the class BDI. The Kitaev honeycomb model [42] can be mapped onto free spinless fermions withp-wave pairing on a hon- eycomb lattice, using the Jordan-Wigner transformation. The solutions for the edge modes for a semi- infinite lattice¹ have been studied earlier [49–52]. The momentum space Hamiltonian in terms of the

1We would like to point out that this system is different from two-dimensionalpx+ipyfermionic superfluids, whose excitation spectra include gapless Majorana-Weyl fermions [43]. Volovik showed that in such chiral superfluids, the fermionic zero modes along the domain wall have the same origin as the fermion zero modes appearing in the spectrum of the Caroli-de Gennes-Matricon bound states in a vortex core [44–48]. This correspondence can be understood by picturing the chiral fermions as orbiting around the vortex axis, analogous to the motion along a closed domain boundary.

Majorana operators is H_h=X

k

³ ˆ a^†

k bˆ^†

k

´ h_h(_k)

µ aˆ_k bˆ_k

¶

, h_h(_k)=

µ 0 A(_k) B(_k) 0

¶

, _k=(kx,ky) ,

A(_k)= −2i

·

J₃+J₁cos

µkx−ky

2

¶ +J₂cos

µkx+ky

2

¶¸

+2

· J₁sin

µkx−ky

2

¶ +J₂sin

µkx+ky

2

¶¸

,

B(k)=2i

·

J3+J1cos

µk_x−k_y 2

¶ +J2cos

µk_x+k_y 2

¶¸

+2

· J1sin

µk_x−k_y 2

¶ +J2sin

µk_x+k_y 2

¶¸

(4.31) with the eigenvalues

E(_k)= ±2 (·

J₃+J₁cos

µkx−ky

2

¶ +J₂cos

µkx+ky

2

¶¸2

+

· J₁sin

µkx−ky

2

¶ +J₂sin

µkx+ky

2

¶¸2)1/2

.^(4.32) We will consider two kinds of edges [49, 50], namely, zigzag and armchair, which can support Majo- rana fermions. We will find the phase diagram using the EP’s corresponding to these edges setting either A(_k)=0^orB(_k)=0, after complexifying the momentum component perpendicular to the edge. For this 2d case,f in equation (2.1) is a function of(J1,J2,J3,k_∥)^{, where}k_∥is the momentum along the 1d edge being considered.

One can have a zigzag edge in they-direction, according to the convention of Nakada et al. [49], so that we will complexifyk_⊥=k_x^{, and}k_∥=k_y will be one of the parameters determining the topological phase transition points. The solution forB(k_x=k_⊥,k_y=k_∥)=0is given by

k_⊥= −2 i ln

·

−J1exp(ik_∥/2)+J2exp(−ik_∥/2) J3

¸

. ^(4.33)

(a) (b)

(c) (d)

Figure 5. (Color online) Topological phase diagrams for edges for the 2d honeycomb lattice described by equation (4.31), as captured by the functionf(J₁,J₂,J₃,k_∥)defined in equation (2.1). Panels (a) and (b) show the number of chiral Majorana zero modes for a zigzag edge, while panels (c) and (d) show the same for an armchair edge located at the top of a semi-infinite lattice.

Majorana zero modes exist for all the values ofk_y^ifJ₁+J₂<J₃. There is no edge state if|J₁−J₂| >J₃^{. For} J₁=J₂=J₃, edge states exist if|k_y| >2π/3. These results have been plotted in figures 5 (a) and 5 (b).

For the armchair edge [49] in thex-direction on the top of the lattice, we will complexifyk_⊥=k_y^{, and} k_∥=k_x will be now one of the parameters determining the topological phase transition points. The two EP’s forA(k_x=k_∥,k_y=k_⊥)=0are given by

k_⊥^±= −2 i ln

·−J₃exp(−ik_∥/2)± q

J₃²exp(−ik_∥)−4J₁J₂ 2J₂

¸

. ^(4.34)

No Majorana zero mode exists for any value ofk_x ^if J₁=J₂=J₃^or J₁<J₂^{. For} J₁>J₂, a Majorana fermion can exist for a specific range of values fork_x. Figures 5 (c) and 5 (d) show these topological phases, obtained using equation (2.1).

Equations (4.33) and (4.34) are seen to coincide with the solutions of the Majorana edge states ob- tained earlier by the transfer matrix formalism [49, 52].

5. EP formalism for the DIII class

A point defect in class DIII can support a Majorana Kramers pair (MKP) corresponding to doubly de- generate Majorana zero modes, whereas a line defect can support a pair of helical Majorana edge states.

Both are characterized by aZ2topological invariant. The chiral symmetry operatorO can be defined such that the Hamiltonian in class DIII can be brought to the block off-diagonal form [equation (2.2)], just like for the class BDI.

5.1. 1d model

A simple 1d model of topological superconductivity in the class DIII is described by the Hamiltonian [53]

H_m1=X

k

Ψ^†_kh_m1(k)Ψk, Ψk=(c_k↑,c_k↓,c^†

−k↓,−c^†

−k↑)^T,

h_m1(k)=[ξm1(k)σ0+λ^Rsin(k)σz]τz+∆cos(k)σ0τx, ξm1(k)=tcos(k)−µ. ^(5.1) This system may be realized in a Rashba wire that is proximity-coupled to a nodelesss_±wave supercon- ductor. The energy eigenvalues are given by:

E₁(k)= ± q

[ξm1(k)−λ^Rsin(k) ]²+∆²cos²(k) , E₂(k)= ± q

[ξm1(k)+λ^Rsin(k) ]²+∆²cos²(k) , ^(5.2) whose plots are shown in figure 6 (a) for∆=0.1t^andλ^R=2t^{, as}µ/tis varied along the horizontal axis.

Observing that a chiral symmetry operatorO =σ0τy exists in the presence ofMz, we rotate the Hamiltonian in equation (5.1) to the chiral basis, where it takes the form

H_m1^chi(k)=U₄^†h_m1(k)U₄=

µ 0 h^u_m1(k) h_m1^l (k) 0

¶

, U₄= 1

p2









0 −i 0 i

0 1 0 1

−i 0 i 0

1 0 1 0







 ,

h^u_m1(k)=^diag

³

i∆cos(k)−ξm1(k)+λ^Rsin(k), i∆cos(k)−ξm1(k)−λ^Rsin(k)

´ , h_m1^l (k)=^diag

³

−i∆cos(k)−ξm1(k)+λ^Rsin(k), −i∆cos(k)−ξm1(k)−λ^Rsin(k)´

. ^(5.3)

The solutions for the EP’s are then given by either det£

h^u_m1(k)¤

=0 ⇒ i∆cos(k)−ξm1(k)=s₁λ^Rsin(k)

⇒ k=k^m1u_s1,s2= −i ln

·µ+s₂q

µ²−λ²R−(t−i∆)² t−i (∆+s₁λ^R)

¸

, ^(5.4)