Quasiclassical theory of disordered multi-channel Majorana quantum wires
QQuasiclassical theory of disordered multi-channel Majorana quantum wires
Patrick Neven, Dmitry Bagrets and Alexander Altland
Institut f¨ur Theoretische Physik, Universit¨at zu K¨oln, K¨oln, 50937, Germany (Dated: October 30, 2018)Multi-channel spin-orbit quantum wires, when subjected to a magnetic field and proximity cou-pled to s -wave superconductor, may support Majorana states. We study what happens to thesesystems in the presence of disorder. Inspired by the widely established theoretical methods ofmesoscopic superconductivity, we develop ´a la Eilenberger a quasiclassical approach to topologicalnanowires valid in the limit of strong spin-orbit coupling. We find that the “Majorana number” M ,distinguishing between the state with Majorana fermion (symmetry class B) and no Majorana (classD), is given by the product of two Pfaffians of gapped quasiclassical Green’s functions fixed by rightand left terminals connected to the wire. A numerical solution of the Eilenberger equations revealsthat the class D disordered quantum wires are prone to the formation of the zero-energy anomaly(class D impurity spectral peak) in the local density of states which shares the key features of theMajorana peak. In this way we confirm the robustness of our previous conclusions [Phys. Rev. Lett. , 227005 (2012)] on a more restrictive system setup. Generally speaking, we find that the qua-siclassical approach provides a highly efficient means to address disordered class D superconductorsboth in the presence and absence of topological structures. PACS numbers: 74.78.Na, 71.23.-k, 73.63.Nm, 74.45.+c
I. INTRODUCTION
Semiconductor quantum wires proximity coupled toa conventional superconductor and subject to a mag-netic field may support Majorana fermion edge states .Building on relatively conventional device technology,this proposed realization of the otherwise evasive Ma-jorana fermion has triggered a wave of theoretical andexperimental activity, which culminated in the recent re-port of a successful observation by several experimentalgroups . In these experiments, evidence for the pres-ence of a Majorana is drawn from the observation of azero bias peak in the tunneling conductance into the wire.While the observed signal appears to be naturally ex-plained in terms of a Majorana end state, two of us havepointed out that a midgap peak might be generated byan unrelated mechanism : in the presence of even verymoderate amounts of disorder, the semiconductor wiresupports a zero energy “spectral peak” (an accumulationof spectral weight at zero energy) which resembles theMajorana peak in practically all relevant aspects. Specif-ically, it is (i) rigidly locked to zero energy, (ii) is of nar-row width of O ( δ ), where δ is the single particle levelspacing, (iii) carries integrated spectral weight O (1), and(iv) relies on parametric conditions (with regard to spinorbit interaction, proximity coupling, magnetic field, andchemical potential) identical to those required by Majo-rana state formation. What makes the spectral peak dis-tinct from the Majorana is that it relies on the presenceof a moderate amount of disorder, viz. impurity scatter-ing strong enough to couple neighboring single particleAndreev levels. Besides, the spectral peak is vulnerableto temperature induced dephasing. While this marks adifference to the robust Majorana state, the reported ex-perimental data does show strong sensitivity to temper-ature, which may either be due to an intrinsic sensitivity of the peak, or due to a temperature induced diminishingof the measurement sensitivity, or both. In either case,the situation looks inconclusive in this regard.Generally speaking, the results of Ref. 6, as well asthose of Refs. 7 and 8 suggest that the observation of amidgap anomaly in the tunneling conductance might bedue to either mechanism, disorder peak, Majorana peak,or a superposition of the two, and this calls for furtherresearch.Our previous study was based on an analyticallytractable idealization of a semiconductor quantum wiresubject to a magnetic field sweep. In the present paper,we will explore the role of disorder scattering within amodel much closer to the experimental setup. The priceto be payed for this more realistic description is that afully analytic treatment is out of the question. Instead,we will employ a semi-analytic approach based on theformalism of quasiclassical Green functions. Introducedin the late sixties , the latter has become an indispens-able tool in mesoscopic superconductivity and quan-tum transport in general . We here argue that quasiclas-sical methods are, in fact, tailor made to the modeling ofMajorana quantum wires (or, more generally, quasi one-dimensional topological superconductors.) Specifically,we will show that in the problem at hand, the quasiclas-sical “approximation” is actually very mild. Further, thequasiclassical Green’s function affords a convenient de-scription of the topological signatures of the system interms of Pfaffians. Finally, the approach can be appliedto systems for a given realization of the disorder, and atnumerical cost much lower than that of exact diagonal-ization approaches. As a result, we will be in a positionto accurately describe local spectral properties within areasonably realistic model of a topological multi-channelsuperconductor. As we are going to discuss below, ourfindings support the principal statement made in Ref. 6. a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b The rest of the paper is organized as follows. In sec-tion II, we discuss the principal role played by disorderin the system, the idea of the quasiclassical approach,and its main results. In section III we specify our modelsystem, the quasi-classical approach is introduced in sec-tion IV, and in V we discuss the numerical solution ofthe quasiclassical equations. We conclude in section VI.A number of technical details are relegated to Appen-dices.
II. QUALTIATIVE DISCUSSION AND RESULTS
A schematic of the device currently under experimen-tal investigation is shown in Fig. 1. A semiconductorquantum wire subjected to strong spin orbit interactionis brought in contact to a superconductor (S), and, via atunnel barrier (T) to a normal metal lead (N). The appli-cation of a small excess voltage, V , to the latter inducesa tunnel current into the central region. The differen-tial conductance dI/dV probes the (tunneling) densityof states at an energy V (units e = (cid:126) = c = 1 through-out) relative to the systems chemical potential, µ . Thephysics we are interested in is contained in a band center( V = 0) anomaly in that quantity.In a manner to be discussed in more detail below, onecontribution to the band center density of states is pro-vided by a Majorana bound state localized at the tunnelbarrier. The second, spectral contribution is generatedby a conspiracy all other low lying quasiparticle states inthe system. Technically, these are Andreev bound statesforming at energies ± E j in the region between the tun-nel barrier and the superconductor. The number of thesestates increases with the extension of the wire – a fewhundred nanometers in the experiment – and the numberof transverse channels below the chemical potential. Inthe presence of disorder, the entity of these states definesan effective “quantum dot”. The proximity of a super-conductor, and the breaking of both spin rotation andtime reversal symmetry imply that the system belongsto symmetry class D .The symmetry of class D random systems implies aclustering of levels at zero energy. Loosely speaking, theconventional level repulsion of random spectra turns intoa zero energy level attraction. On a resolution limitedto scales of order of the mean level spacing, the zeroenergy density of states (DoS) is enhanced by a factorof two relative to the mean background, i.e. it showsa peak. This phenomenon manifests itself at relativelysmall sample-to-sample fluctuations, i.e. the peak is asample specific effect. In passing we note that the weakanti-localization phenomenon discussed in Ref. 8 restson the same principal mechanism of midgap quantuminterference.Below, we will explore the phenomenon of spectralpeak formation in a setting modelled to closely mimic the“experimental reality”. To be more specific, we considera semiconductor wire supporting a number of N > xzN T SL ∆ VµB
FIG. 1. Disordered spin-orbit nanowire, subjected to colinearmagnetic field, is proximity coupled to the s -wave supercon-ductor (S) and terminated by the tunneling barrier (T) at oneof its ends. A sketch below shows the profile of the inducedsuperconducting gap ∆( x ) and gate induced potential V ( x )defining the tunnel barrier. Andreev bound states (ABS) aredepicted by dotted lines. transverse channels below the chemical potential. Weassume a value of the chemical potential such that thehighest lying of these channels is “topological” (chemicalpotential falling into the gap opened by the simultaneouspresence of order parameter and magnetic field.)To quantitatively describe this phenomenon, we gen-eralize the quasiclassical Green function approach to thepresent setting. Indeed, the phenomena we are interestedin manifest themselves on length scales large in compar-ison to the Fermi wavelength, yet smaller than the rele-vant coherence length, and this makes them suitable toquasiclassical treatment. The quasiclassical theory of thepresent paper is formulated in terms of the Eilenbergerfunction ˜ Q ( x ; (cid:15) ), which is a position and energy depen-dent matrix of size 8 N × N . The factor of 8 = 2 × × Q ) will be seen to assume two values( ± Z invariantof the underlying one-dimensional class D superconduc-tor, in dependence on system parameters: for values ofthe magnetic field smaller or larger than a critical field, B c ≡ (cid:112) ∆ + µ , where µ and ∆ are chemical poten-tial and bulk order parameter, the invariant is trivial ornon-trivial, resp. In the latter case, the wire supports aMajorana state at the barrier, in the former it doesn’t.Fig. 2 shows a profile of the local DoS (LDoS) com-puted at the left end of the wire for typical system pa-rameters as detailed below. The green dashed and redsolid curve correspond to a situation without and withMajorana state. Within our model, the states in the wirecarry a narrow width, ∼ g T δ , reflecting the possibility ofdecay through the tunnel barrier into the adjacent lead. μ = 0.5Δ, B = 0.25Δμ = Δ, B = 2.66Δ < D o S > [ ν ] ε /Δ −1 −0,5 0,5 1 ν ε FIG. 2. Disorder averaged local density of states (LDoS) atthe left end of the spin-orbit quantum wire sketched in Fig. 1(in units ν = 1 / πv ). The number of occupied bands N = 2corresponds to four transport channels. Parameters are: (i)red solid line ( M = − B = 2 . µ = ∆ — topologicalphase; (ii) green dashed line ( M = +1), B = 0 . µ = 0 . L = 4 v/ ∆, dimensionlessstrength of disorder γ w /v = 0 . l = 0 . L . Tunneling rate Γ = 5 · − ∆.Velocities in two bands were taken to be equal, v = v = v .The inset shows profiles of the DoS resulting from randommatrix theory. (Here, δ is the mean level spacing of the wire, and g T (cid:46) negative spectral peak atzero energy (the negative of the positive class D peak),superimposed on a single δ -peak (the Majorana). Thejoint signature of these two structures is seen in the solidcurve in Fig. 2.At resolutions limited to values ∼ δ , the superpositionof the Majorana and the class B peak looks next to in-distinguishable from the class D peak (dashed), and thissimilarity of unrelated structures might interfere with theunambiguous observation of the Majorana by tunnelingspectroscopy. Indeed, the differential tunneling conduc-tance monitored in experiment, dIdV = e g T π (cid:126) (cid:90) ∞−∞ ∂f F ∂(cid:15) ( (cid:15) − V ) ν L ( (cid:15) ) ν d (cid:15), (1) D o S [ ν ] ε /Δ −1 −0,5 0,5 1 ν ε FIG. 3. The average LDoS (solid green line) and the squareroot of its second (reducible) moment (cid:104) ν L ( (cid:15) ) (cid:105) / (dashed blueline) at the left end of the spin-orbit quantum wire in thetrivial phase ( M = +1). System parameters are listed inFig. 2. The inset shows profiles of the mean DoS and thesquare root of the two level correlation function resulting fromrandom matrix theory. is essentially determined by the local density of states,i.e. the structures shown in Fig. 2 are expected to reflectdirectly in the measured signal. (Here, f F is the Fermidistribution, ν L is the DoS at the left barrier, ν is theDoS in the single chiral channel per unit length, and V the applied voltage.)The profiles of the curves shown in Fig. 2 were com-puted for a two-channel quantum wire at a mean freepath l (cid:39) L of the order of the system size. We are ad-dressing a system at the interface between the ballisticand the localized regime. In view of these system param-eters it is remarkable that the DoS profiles in Fig. 2 showstriking similarity to the average DoS of a class D andB random matrix model . For comparison the averageDoS of a class B and D random matrix Hamiltonian isshown as an inset in Fig. 2. The similarity of the resultsindicates that the system of subgap states in our systembehaves as if it formed an effective chaotic quantum dotlocalized in the vicinity of the left system boundary (rightto the tunnel barrier). Taking into account spin, chiraland channel quantum numbers, the mean level spacingin such dot is given by δ (cid:39) πv/ N L . In the presenceof magnetic field the BCS gap in the superconducting re-gion of the wire reads (cid:15) − = | B − (cid:112) ∆ + µ | . The numberof subgap Andreev levels forming the effective dot is thusgiven by N levels (cid:39) (cid:15) − /δ .The profiles shown in Fig. 2 are ensemble averages, (cid:104) ν L (cid:105) of the LDoS, ν L , where the sampling was over ∼ D o S [ ν ] ε /Δ −1 −0,5 0 0,5 1 FIG. 4. The sample specific LDoS in the trivial phase withoutMajorana state ( M = +1) for two different disorder realiza-tions. Tunneling rate Γ = 0 . ± (cid:15) min lyingclose to Fermi energy and having energy splitting ∼ Γ (dottedgreen line); (ii) particle and hole states have merged into asingle zero-energy peak of width Γ and can not be resolvedby tunnel spectroscopy (solid magenta line). For the chosenset of parameters the mean level spacing δ = 0 .
2∆ and thegap in the S region (cid:15) − = 0 . N levels (cid:39) average (cid:112) (cid:104) ν L (cid:105) . The relatively minor deviation betweenaverage and typical DoS demonstrates that the standarddeviation δν ( (cid:15) ) = (cid:68)(cid:16) ν ( (cid:15) ) − (cid:104) ν ( (cid:15) ) (cid:105) (cid:17) (cid:69) / (2)characterizing the strength of mesoscopic fluctuations isrelatively small.The weakness of fluctuations implies that the disor-der peak is a realization specific phenomenon. Thisis demonstrated explicitly in Fig. 4, where two un-averaged DoS profiles individual disorder configurationsare shown. The data is for a “non-topological” system,no non-degenerate zero energy level is present. Depend-ing on whether the pair of lowest lying Andreev states(+ (cid:15) min , − (cid:15) min ) exceeds the level broadening Γ, the DoSenhancement will assume the form of a single peak (ma-genta solid), or a split peak (green dashed). In eithercase, an excess DoS of integrated spectral weight ∼ III. MODEL OF THE MAJORANA NANOWIRE
We consider a multi-band quantum wire of width L z ,subject to Rashba spin-orbit coupling, proximity cou-pling to an s -wave superconductor, and to a magneticfield . We choose coordinates such that the wire liesalong the x -axis, parallel to the magnetic field, the y -axisis perpendicular to the surface of the superconductor, andthe spin orbit field is pointing along the z -axis.In the rest of this section we will specify theBogoliubov-de-Gennes (BdG) Hamiltonian describingthis system in the presence of disorder (section III A).We will then linearize the electron spectrum in the limitof strong spin-orbit coupling thereby introducing a de-scription in terms of right and left one-dimensional chiralfermions (section III B), and finally transform the Hamil-tonian into Majorana basis (section III C). Experts in thedescription of topological quantum wires may proceed di-rectly to the end of the section. A. Bogoliubov-de-Gennes Hamiltonian
Introducing a four component spinor Ψ =( ψ ↑ , ψ ↓ , ¯ ψ ↑ , ¯ ψ ↓ ) in the product of spin and particle-hole spaces, the BdG Hamiltonian H describing thesystem in the xz -plane reads H = 12 (cid:90) ¯Ψ( x, z ) (cid:18) ˆ h + ˆ W i ˆ s y ∆ ∗ − i ˆ s y ∆ − ˆ h T0 − ˆ W T (cid:19) Ψ( x, z )d x d z, (3)withˆ h = − ( ∂ x + ∂ z ) / m − µ ( x ) + B ( x )ˆ s x − iα (ˆ s z ∂ x − ˆ s x ∂ z ) . Here, ∆ = ∆( x ) is the proximity amplitude induced thesuperconductor, B ≡ gµ B H , where H is the externalfield, we have taken into account the transverse momen-tum ( − i∂ z ) in the Rashba term, and ˆ W ( x ) is the randomdisorder Hamiltonian. The Pauli matrices ˆ s operate inspin space.We proceed by introducing a system of transverse wavefunctions, { Φ σn ( z ) } , which leads to a linear decompositionof the Grassmann fields ( σ = ↑ , ↓ ) as ψ σ = (cid:88) nσ Φ σn ( z ) ψ ( n ) σ ( x ) , ¯ ψ σ = (cid:88) nσ Φ σn ( z ) ∗ ¯ ψ ( n ) σ ( x ) . (4)Defining Ψ ( n ) = ( ψ ( n ) ↑ , ψ ( n ) ↓ , ¯ ψ ( n ) ↑ , ¯ ψ ( n ) ↓ ) T , ¯Ψ ( n ) = ( ¯ ψ ( n ) ↑ , ¯ ψ ( n ) ↓ , ψ ( n ) ↑ , ψ ( n ) ↓ ) , (5)the Hamiltonian then takes the form H = 12 (cid:90) ¯Ψ ( n ) ( x ) (cid:32) ˆ h ( n )0 δ nm + ( iα ˆ s x ∂ z ) nm + ˆ W nm i ˆ s y ∆ ∗ δ nm − i ˆ s y ∆ δ nm − (ˆ h ( n )0 ) T δ nm − ( iα ˆ s x ∂ z ) mn − ˆ W T mn (cid:33) Ψ ( m ) ( x )d x, (6)where the one-dimensional Hamitonian ˆ h ( n )0 acts in the n -th band asˆ h ( n )0 = − ∂ x / m + µ ( n ) z − µ ( x ) + B ( x )ˆ s x − iα ˆ s z ∂ x . (7)In the case of an ideal waveguide the transverse wavefunc-tions are Φ σn ( z ) = (cid:112) /L z sin( nπz/L z ), so that µ ( n ) z = π n / (2 mL z ) and the matrix elements of the spin-orbitinteraction read h s . o .nm = ( iα ˆ s x ∂ z ) nm = − iαL z nmn − m (1 − ( − n + m ) ˆ s x . (8)Let us now assume a thin wire, L z (cid:46) l so = (cid:126) / ( mα ),where l so is the spin-orbit length. In this case the ma- trix elements h s . o .nm (cid:28) µ ( n ) z can be treated as perturba-tions. To this end we introduce a unitary transforma-tion ˆ U , which brings the high energy part of the Hamil-tonian, ( µ ( n ) z δ nm + h ( n ) nm ), to diagonal form. Due tothe weakness of the perturbation, the transformation isclose to unity, ˆ U = exp( i ˆ X ) (cid:39) i ˆ X , with generatorsˆ X nm = O ( L z /l so ). To first-order perturbation theorytheir explicit form readsˆ X n (cid:54) = m (cid:39) ih s . o .nm µ ( m ) z − µ ( n ) z , ˆ X nn = 0 . (9)The transformation Ψ → ˆ U Ψ and ¯Ψ → ¯Ψ ˆ U † generatesthe approximate Hamiltonian H (cid:39) (cid:90) ¯Ψ ( n ) ( x ) (cid:32) ˆ h ( n )0 δ nm + ˆ V nm + ˆ W nm i ˆ s y ∆ δ nm + ( δ ˆ∆) nm − i ˆ s y ∆ δ nm − ( δ ˆ∆) ∗ nm − (ˆ h ( n )0 ) T δ nm − ˆ V T mn − ˆ W T mn (cid:33) Ψ ( m ) ( x ) d x. (10)Here we have redefined our notation for the disorderpotential ˆ W → ˆ U ˆ W ˆ U † , and introduced a small correc-tion to the quasiparticle Hamiltonian, ˆ V = α [ ˆ X, ˆ s z ] ∂ x ,and to the order parameter, δ ˆ∆ = − ∆[ ˆ X, ˆ s y ]. Onecan estimate the matrix elements of these operators as V nm ∼ (cid:15) ( L z /l so ) and ( δ ∆) nm ∼ ∆( L z /l so ), where (cid:15) isthe quasiparticle energy. Since we are primary interestedin the effects of disorder, we will limit our discussionto the situation when the low-energy physics is disorderdominated, i.e. the random potential ˆ W > max { ˆ V , δ ˆ∆ } masks the off-diagonal matrix elements of the determinis-tic Hamiltonian. Having this in mind we thus omit bothˆ V and δ ˆ∆ throughout the paper. However, this assump-tion does not affect our results in qualitative ways.Assume that the chemical potential lies close to thebottom of the N -th band, µ (cid:39) µ ( N ) z . We refer to thisband as the “topological band”, since in the absence ofinterband scattering it defines whether the quantum wireis in the topological phase or not. The condition of thetopologically non-trivial phase reads B > B c = ∆ +( µ − µ ( N ) z ) . We also assume a hierarchy of energy scales µ ( N ) z − µ ( N − z (cid:29) E so (cid:29) ∆ ∼ B , where E so = mα / n < N are in thetrivial phase, since for them B (cid:28) ∆ + ( µ − µ ( n ) z ) . Inthe current experiments E so (cid:38) ∆, i.e. our theory is onthe border of applicability (see, however, the discussionin section IV). B. One-dimensional chiral fermions
For strong spin-orbit coupling, E so (cid:29) B , each band ischaracterized by two Fermi momenta, k a/bn = ∓ mα + (cid:114) m (cid:16) µ ( N ) z − µ ( n ) z + E so (cid:17) , (11)as shown in Fig. 5. We aim to construct a low energyHamiltonian describing the system at energy scales (cid:46) E so . To this end we represent the fields ψ ( n ) σ as ψ ( n ) ↑ ( x ) (cid:39) R ( n ) a ↑ ( x ) e ik an x + L ( n ) b ↑ ( x ) e − ik bn x , (12) ψ ( n ) ↓ ( x ) (cid:39) L ( n ) a ↓ ( x ) e − ik an x + R ( n ) b ↓ ( x ) e ik bn x , (13)in terms of a superposition of right ( R ) and left ( L )chiral fermions, and then linearize the Hamiltonian (10) FIG. 5. Chiral one-dimensional fermions in the spin-orbitquantum wire: a and b channels are shown by open and filleddots, resp. The magnetic field B affects the dispersion relationfor the a -channel only if the chemical potential is close to µ ( N ) z . around the Fermi momenta. Modes with Fermi momentaof equal modulus define a conduction channel. We ob-serve that each channel belongs to one of the two subsets, a or b . In the a -channel R -movers have spin up, and L -movers have spin down, while for the channels of type b spins are reversed. The a -channel of band index N playsa distinguished role in that it has zero Fermi momentum, k aN = 0. This particular mode defines what we call the“topological channel”. For max { B, ∆ } (cid:28) E so it is theonly channel strongly susceptible to the magnetic field,thus defining the topological phase of the whole wire. The effect of B (but not ∆) on other channels can besafely neglected.We next collect all R and L fields into two spinors,Ψ = ( R a ↑ , R b ↓ , L a ↓ , L b ↑ , ¯ R a ↑ , ¯ R b ↓ , ¯ L a ↓ , ¯ L b ↑ ) , (14)and ¯Ψ = ( σ ph x Ψ) T where the band index ( n ) is left im-plicit. The spinor Ψ acts in an 8 × N -dimensional space,defined by the direct product of band (index n ), chan-nel ( a/b ), chiral ( R/L ) and particle-hole spaces. In thisrepresentation the low energy Hamiltonian is given by H = N (cid:88) n,m =1 (cid:90) ¯Ψ ( n ) ( x ) (cid:32) ˆ h ( n )0 δ nm + ˆ W nm i ∆ σ RLy ⊗ σ abz δ nm − i ∆ σ RLy ⊗ σ abz δ nm − (ˆ h ( n )0 ) T δ nm − ˆ W T mn (cid:33) ph Ψ ( m ) ( x ) d x, (15)where ˆ h ( n )0 = (cid:18) − iv n ∂ x − µ ( B/
2) ( ab + σ abz ) δ Nn ( B/
2) ( ab + σ abz ) δ Nn iv n ∂ x − µ (cid:19) RL . (16)Here, v n = [2 m ( µ ( N ) z − µ ( n ) z + E so )] / /m is the Fermivelocity of the n -th channel, the chemical potential µ isnow defined relative to the energy µ ( N ) z , and ˆ W nm ( x ) isa random 8 × W nm ) † = ˆ W mn . In deriving this Hamiltonian wehave neglected oscillatory terms with phases ix ( ± k a ± k b ).This is justified because the typical lengths involved inthe problem ( ξ ∼ ∆ /v and l B ∼ B/α ) exceed by far theFermi length λ F ∼ max(1 /k a , /k b ) determined by thelarge spin-orbit energy E so . We have also used that theorder parameter ∆( x ) can be chosen to be real. C. Majorana representation
The 8 N × N first-quantized matrix defining theHamiltonian (15) satisfies to the p-h symmetry,ˆ H = − σ ph x ˆ H T σ ph x , (17)which is the defining condition for a class D Hamilto-nian (here, the transposition acts on kinetic term as ∂ T x = − ∂ x .) However, all what follows will be more con- veniently formulated in an alternative representation, inwhich the symmetry assumes a different form: for eachband n we define a set of eight Majorana fields ξ Ra = ( R a ↑ + ¯ R a ↑ ) / √ , η Ra = ( R a ↑ − ¯ R a ↑ ) / √ i,ξ Rb = ( R b ↓ + ¯ R b ↓ ) / √ , η Rb = ( R b ↓ − ¯ R b ↓ ) / √ i, (and analogous relations for the L –movers, with spin re-versed), which we combine into the 8 N –spinor˜ χ = ( ξ Ra , ξ Rb , ξ La , ξ Lb , η Ra , η Rb , η La , η Lb ) T . (18)The two spinors χ and Ψ are related by the unitary trans-formation ˜ χ = U Ψ , ˜ χ T = ¯Ψ U † , (19) U = ab ⊗ RL ⊗ √ (cid:18) i − i (cid:19) ph . We combine the Majorana fields of a and b channels as ξ R = ( ξ Ra , ξ Rb ) T , and reorder the spinor components as, χ = ( ξ R , η L , η R , ξ L ) T . (20)In this representation (which will be used throughout therest of the paper) the Hamiltonian (15) takes the form H = N (cid:88) n,m =1 (cid:90) χ T n ( x ) ˜ H nm ( x ) χ m ( x ) d x, ˜ H nm ( x ) ≡ (cid:32) ˆ h ( n ) − δ nm + i ˆ W −− nm iµ σ RLz ⊗ ab δ nm + i ˆ W − + nm − iµ σ RLz ⊗ ab δ nm + i ˆ W + − nm ˆ h ( n )+ δ nm + i ˆ W ++ nm (cid:33) , (21)where the deterministic part readsˆ h ( n ) ± = − iv n σ RLz ⊗ ab ∂ x − σ RLy ⊗ ˆ∆ ( n ) ± , (22)∆ ( n ) ± = ∆ σ abz ± ( B/ ab + σ abz ) δ Nn , (23) and the random matrices are constructed asˆ W −− = (cid:18) ˆ w RR − ˆ w RL ˆ w LR ˆ w LL (cid:19) ab , ˆ W − + = (cid:18) − ˆ w RR ˆ w RL ˆ w LR ˆ w LL (cid:19) ab , (24)in terms of real (symmetric) and imaginary (antisymmet-ric) parts of the random matrices ( ˆ W ) RL = ˆ w RL + i ˆ w RL etc. The remaining blocks of the ˆ W -matrix are definedby ˆ W + − = − σ abz ˆ W − + σ abz , ˆ W ++ = σ abz ˆ W −− σ abz . (25)In the Majorana basis (20), the class D symmetry is ex-pressed through the antisymmetry ˆ H T = − ˆ H .We finally specify the statistics of disorder. We chooseˆ W ( x ) to be a δ -correlated and Gaussian distributed ran-dom matrix of size 8 N × N with a zero mean value (cid:104) ˆ W ( x ) (cid:105) and variance (cid:104) w ij ( x ) w i (cid:48) j (cid:48) ( x (cid:48) ) (cid:105) = γ w δ ( x − x (cid:48) )( δ ii (cid:48) δ j (cid:48) j (cid:48) + δ ij (cid:48) δ ji (cid:48) ) , (cid:104) w ij ( x ) w i (cid:48) j (cid:48) ( x (cid:48) ) (cid:105) = γ w δ ( x − x (cid:48) )( δ ii (cid:48) δ jj (cid:48) − δ ij (cid:48) δ ji (cid:48) ) . (26)Here, the composite indices ( i, j etc.) label states inthe direct product of band, channel and chiral spaces.The scattering matrices defined in this way break time-reversal and spin rotation symmetry (e.g., random spin-flip scattering caused by random spin-orbit terms is in-cluded in (26).) The strength of disorder set by the co-efficient γ w translates into the “golden rule” scatteringrate τ − = 2 γ w N (cid:88) n =1 (1 /v n ) (27)of the normal conducting (i.e. superconductor decou-pled) quantum wire. IV. QUASICLASSICAL APPROACH
The kinetic term in the low energy Hamiltonian (21)which was discussed in the previous section is linear inmomentum, and this facilitates the formulation of quasi-classical equations of motion (aka Eilenberger equations)for the model at hand . We here review the construc-tion of these equations in a manner closely following thespirit of Refs. 19 and 20. After introducing the basicsof the method (section IV A) we construct the Eilen-berger Q –function in the limit of a single clean “topo-logical channel” (section IV B) and discuss the resultingdensity of states (section IV C). In section IV D we definethe Z topological invariant in terms of the Q –matrix.Section IV E outlines the general construction of the so-lution to the Eilenberger equation in the inhomogeneousdisordered wire with boundary conditions. A. Eilenberger method
We start by defining, g R/A(cid:15),nm ( x, x (cid:48) ) = √ v n G R/A(cid:15),nm ( x, x (cid:48) ) σ RLz √ v m , (28) where G R/A(cid:15) ( x, x (cid:48) ) ≡ (cid:104) x | ( (cid:15) ± i − ˜ H ) − | x (cid:48) (cid:105) are the retardedand advanced Green’s functions of the system. Undertransposition (which in our current representation repre-sents the particle-hole symmetry) the function g behavesas g R(cid:15) ( x, x (cid:48) ) = − σ RLz (cid:2) g A − (cid:15) ( x (cid:48) , x ) (cid:3) T σ RLz , (29)i.e. advanced and retarded Green’s functions get inter-changed. It is not hard to derive the two mutually adjointdifferential equations ∂ x g (cid:15) ( x, x (cid:48) ) + L (cid:15) g (cid:15) ( x, x (cid:48) ) = − iδ ( x − x (cid:48) ) ,∂ x (cid:48) g (cid:15) ( x, x (cid:48) ) − g (cid:15) ( x, x (cid:48) ) L (cid:15) = iδ ( x − x (cid:48) ) , (30)describing the dynamical evolution of g . Here, L (cid:15) ≡ − i (ˆ ω − P ) (31)where matrix ˆ ω has elements(ˆ ω ) nm ≡ ( (cid:15)/v n ) σ RLz δ nm (32)and the operator P is related to the Hamiltonian matrix˜ H as P nm = i∂ x + v − / n ( σ RL ˜ H nm ) v − / m . (33)Due to the antisymmetry of ˜ H = − ˜ H T , the operator P obeys the particle-hole symmetry P = − σ RLz P T σ RLz . (34)We next define the Eilenberger function as Q (cid:15) ( x ) = lim x (cid:48) → x (cid:2) i g (cid:15) ( x, x (cid:48) ) − sgn( x − x (cid:48) ) (cid:3) , (35)where the subtraction of the sgn-function regularizes adiscontinuity arising in g at x = x (cid:48) due to the combina-tion of linear derivatives and δ -function inhomogeneityin Eqs. (30). Subtracting the two equations in (30), wethen obtain the Eilenberger equation of motion ∂ x Q (cid:15) ( x ) + (cid:2) L (cid:15) , Q (cid:15) ( x ) (cid:3) = 0 . (36)The Eilenberger function Q obeys the particle-hole sym-metry σ RLz Q T − (cid:15) ( x ) σ RLz = − Q (cid:15) ( x ) , (37)and the normalization condition Q (cid:15) ( x ) = , where isthe unit matrix (the latter condition can be checked byverifying that Eq. (36) preserve the normalization Q =const . ). The unit-value of the normalization constant isfixed by the jump height of the sgn-function in (35). Wefinally note that the operator P can be straightforwardlyconstructed from (33), however, for our present purposes,we need not state its explicit form in generality. B. Eilenberger function in the clean limit
As a warmup, we apply the quasiclassical approachto the limit ( ˆ W = 0) of an infinite clean quantum wiresubject to constant B, ∆. In this system all channels aredecoupled, and we may concentrate on the 4 × Q (cid:15) ( B, µ ) describing the “topological” channel a with n = N . (The Eilenberger function of the other channels maybe obtained by setting B = 0, rescaling the velocity andtransforming ∆ → − ∆ in the case of b -type channels.)We start by introducing the 4 × L (cid:15) = − i (cid:18) (cid:15) σ RLz − i ∆ − σ RLx − iµ RL iµ RL (cid:15) σ RLz − i ∆ + σ RLx (cid:19) , (38)as the reduction of the general operator L (cid:15) to a singlechannel. The solution Q (cid:15) is then determined by the rela-tions [ Q (cid:15) , L (cid:15) ] = 0 and Q (cid:15) = . To solve these equations,we assume L (cid:15) to be diagonalized as L (cid:15) = T (ˆ λ ⊗ σ RLz ) T − , ˆ λ ≡ (cid:18) λ + λ − (cid:19) , (39)where the exact form of the (non-unitary) transformationmatrix T will not be needed and λ ± = (cid:112) λ ± λ,λ = B + ∆ − µ − (cid:15) ,λ = 2 (cid:112) B ∆ − (∆ − (cid:15) ) µ (40)are the eigenvalues. The defining equations for Q arethen solved by matrices of the form Q (cid:15) = T Λ T − , (41)where Λ = ( ± , . . . , ±
1) is a diagonal 4 × (cid:15) → (cid:15) ± i (cid:15) ± i σ RLz , which means that the ap-propriate matrix structure of the retarded Green’s func-tion (opposite for advanced) is given byΛ = σ RLz . (42)A more explicit derivation of this structure is detailed inAppendix A. C. Clean density of states
The density of states in the bulk of the topological wireis given by ν ( (cid:15) ) = (2 πv ) − Re tr σ RLz Q (cid:15) . The matrices T diagonalizing Q do not commute with σ RLz , which meansthat a little extra work is required to evaluate the trace.We start from the representation Q (cid:15) = L (cid:15) λ + P + + L (cid:15) λ − P − , (43) where P + and P − are projectors on the space of L (cid:15) -eigenstates with eigenvalues ± λ + and ± λ − , resp.: P + = T diag( , T − ,P − = T diag(0 , ) T − . (44)That this representation faithfully represents the matrix Q (cid:15) is checked by application of (43) in the eigenbasiswhere all matrices assume a diagonal form. It remainsto obtain a representation of P ± which does not makeexplicit reference to the diagonalizing matrices T . To thisend, notice that (eigenrepresentation understood) L (cid:15) =diag( λ , λ − ) = λ + λP + − λP − . This equationcan straightforwardly solved as P ± = 12 (cid:18) ± λ ( L (cid:15) − λ ) (cid:19) . (45)Substituting this expression into (43), we obtain a repre-sentation of Q which makes reference only to the operator(38), and its eigenvalues. Computing the trace, we ob-tain DoS profiles as shown in Fig. 6. Before discussingthe structure of these results, a general remark may be inorder: the DoS of a one dimensional quantum system isdetermined by an interplay of the kinetic energy operator( k ↔ − i∂ x ) and the “potential” ( L ). On the other hand,we computed the DoS from Q as determined by L , andthis matrix seems to be oblivious to the kinetic energy.A closer look, however, shows that information on theband dispersion sneaks in via the nonlinear constraint Q = . Indeed, the conservation of the constraint, andthe unit value of the normalization are consequences ofthe linearity of the derivative operator in (30), which inthis way co-determines the structure of Q .Inspection of Eq. (43) shows that the DoS containssingularities at the zeros λ ± ( (cid:15) ) = 0, which are located at (cid:15) ± = | B ± (cid:112) ∆ + µ | . (46)Let us assume that B > ∆. Fig. 6(a) shows the ensuingDoS profile, along with the underlying dispersion relationfor a value of the chemical potential µ < µ c , where µ c = (cid:112) B − ∆ , (47)defines a critical value where the lower of the DoS sin-gularities, (cid:15) − , touches zero and the band gap closes[Fig. 6(b)]. At larger values µ > µ c , the gap reopens, (c),the DoS looks qualitatively similar to that of the µ < µ c regime, but the system is in a topologically distinct state(see the next section.) Finally, at values µ > µ ∗ , where µ ∗ = (cid:113) B + B (cid:112) B + 4∆ / √ , (48)the lower band (cid:15) − ( k ) develops an extremum at finite val-ues of k which manifests in a third van-Hove singularityat the energy (cid:15) = ∆ (cid:112) − B /µ , (49)as shown in panel (d) of Fig. 6. ν ( (cid:15) ) (cid:15) (a) µ < µ c k(cid:15) k ν ( (cid:15) ) (cid:15) (b) µ = µ c k(cid:15) k ν ( (cid:15) ) (cid:15) (c) µ > µ c k(cid:15) k ν ( (cid:15) ) (cid:15) (d) µ > µ ∗ k(cid:15) k (cid:15) + (cid:15) − (cid:15) (cid:15) + (cid:15) − (cid:15) + (cid:15) − (cid:15) + FIG. 6. Sketch of the density of states and the corresponding dispersion relations in the clean nanowire. The gap is closed atthe transition point µ = µ c (a), and opened again (b,c). For µ > µ ∗ two minima develop (d). D. Z topological invariant The symmetry Eq. (37) implies that at zero energythe product σ RLz Q (cid:15) =0 is an antisymmetric 4 × Q is unity, the sameis true for the determinant of the 4 × σ RLz = σ RLz ⊗ . Consequently, the Pfaffian of σ RLz Q (cid:15) =0 —which squares to the determinant of that matrix — maytake one of two values, ±
1. This motivates the definitionof the topological index N top = Pf( σ RLz Q (cid:15) =0 ) = (cid:26) +1 , µ > µ c − , µ < µ c , (50)distinguishing between the two phases. Computing N top from (43), we find the index structure stated in (50).(Right at the critical point µ = µ c , the matrix Q (cid:15) =0 becomes singular and the index cannot be defined.) E. Eilenberger equation with disorder
In this section we discuss the formal solution of theEilenberger equation in the presence of disorder. The so-lution is “formal” in the sense that the Eilenberger Greenfunction will be a functional of a given realization of thedisorder. To obtain practically useful information, onewill want to average over different realizations, and thisstep of the computation needs to be done numerically, asdiscussed in the next section.To start with, consider the prototypical system geom-etry shown in Fig. 7. The terminals indicated at the leftand right represent superconducting regions, assumednon-disordered for simplicity. (This is an inconsequen-tial assumption provided the rate of disorder scattering τ − ∼ N γ w (cid:104) /v n (cid:105) does not exceed the energy gaps ( (cid:15) − or (cid:15) ) in the terminals.) In these regions, the Eilenbergerequation can be solved analytically, as discussed in theprevious section. Describing the disorder present in thecenter region in terms of a generalized variant of the qua-siclassical evolution operator L , we will show how the leftand right asymptotic configuration of the Green function get connected by a transfer matrix, M , functionally de-pending on the disorder configuration. The ensuing gen-eralized Green function will then be the starting pointfor our numerical analysis.To be more specific, we consider a quantum wire wherethe gap ∆( x ) and/or chemical potential µ ( x ) vary inspace in the region | x | < L/ L/R and µ R/L at x (cid:28) − L/ x (cid:29) L/
2, respec-tively (Fig. 7). These constants set asymptotic values ofthe Q -matrix, Q (cid:15) ( x → −∞ ) ≡ Q − , Q (cid:15) ( x → + ∞ ) ≡ Q + , (51)where Q ± are constructed using the results of sec-tion IV B for the homogeneous profile of ∆, B and µ .The boundary Green’s function Q − and Q + may describedifferent or equivalent topological phases of the wire.We denote the Q -matrices obtained at the interfacebetween the asymptotic superconducting regions, and thecenter disordered region, resp., as Q R = Q ( x R ) , Q L = Q ( x L ) , (52)where x R = − x L = L/
2. These two configurations arerelated by a transfer matrix, Q R = M ( x R , x L ) Q L M − ( x R , x L ) . (53)For arbitrary positions x and x (cid:48) the formal expressionfor the transfer matrix M ( x, x (cid:48) ) at given energy (cid:15) follows xzQ − Q + x L x R Q L Q R FIG. 7. Disordered spin-orbit wire connected to two idealsuperconducting terminals which are described by the Eilen-berger functions Q + and Q − . The Q -matrix at the boundariesbetween the scattering region and terminals is denoted by Q R and Q L . In the superconductors Q –matrix rapidly convergesto either Q + or Q − on a scale of coherence length. M (cid:15) ( x, x (cid:48) ) = P x exp (cid:40) i (cid:90) x (cid:48) x [ˆ ω − P ( y )] d y (cid:41) , (54)with P x denoting the path-ordering operator. A relationsimilar to Eq. (53) connects the boundary matrices Q R/L with the Green’s function in the far right/left region ofthe wire, Q ± = M ( x, x R/L ) Q R/L M − ( x, x R/L ) , (55)assuming that x → ±∞ . The transfer matrix satisfiescertain symmetries. Along with the unitarity condition(cf. Eq. (33)) σ RLz M † (cid:15) σ RLz = M − (cid:15) , (56)it also satisfies to the p-h symmetry: with the use ofEq. (33) one obtains L T (cid:15) = − σ RLz L − (cid:15) σ RLz , and this yields σ RLz M T − (cid:15) σ RLz = M − (cid:15) . (57)We now aim to represent the Q -matrix in the scatteringregion in terms of the transfer matrix M (cid:15) and the asymp-totic Eilenberger functions Q ± . We start by translatingthe transfer matrix relation (55) to a set of algebraic con-ditions relating the matrices Q R/L to Q ± . To this end,notice that the action of the non-unitary (cf. Eq. (56))transfer matrix on a generic matrix Q R will in generalproduce exponentially increasing and decreasing contri-butions. The former are inacceptable in that they leadto exponential divergencies in the quasiclassical Greenfunction. As is detailed in Appendix B, the requirementof a non-divergent Green function leads to the algebraicconditions ( + Q − )( − Q L ) = 0 , (58)( + Q L )( − Q − ) = 0 , while the right matrix Q R should obey the analogousrelations ( − Q + )( + Q R ) = 0 , (59)( − Q R )( + Q + ) = 0 . We finally combined these equations with the transfermatrix relation (55) to obtain closed expressions for Q L/R in terms of the asymptotic configurations Q ± and M .As a result of a straightforward calculation detailed inAppendix B we obtain Q R = + 2 Q + + M Q − M − ( − Q + ) , (60) Q L = + ( − Q − ) 2 Q − + M − Q + M , (61) where M ≡ M ( x R , x L ).These formulae define the starting point for an ef-ficient numerical computation of the disordered Eilen-berger function Q ( x ). To this end, one computes thetransfer matrix M by numerical solution of correspond-ing system of linear first order differential equations. Onenext applies Eqs. (60) and (61) to obtain Q R/L . Finally,our main object of interest, Q ( x ), is obtained by appli-cation of M ( x, x R/L ) to either Q R or Q L .So far we have considered a quantum wire connectedto two superconducting terminals. However, the gener-alization of the method to the system shown in Fig. 1is straightforward. The key observation is that in thelimit of vanishing barrier conductance g T (cid:28) R a ↑ ( x L ) = e iφ L b ↑ ( x L ) , R b ↓ ( x L ) = e iφ L a ↓ ( x L ) , (62)where φ is (energy dependent) reflection phase shift.(Here, we assumed the absence of barrier spin flip scatter-ing, inter-channel scattering or related complications.) Inthe limit of asymptotically high potential barrier φ = π .Relations (62) define boundary conditions for the Eilen-berger function Q L . As verified in Appendix D, theseconditions assume the form of Eqs. (58), where, how-ever, the role of Q − is taken by an “effective” matrix Q − ≡ Q − ( φ ) describing the tunnel junction. The ex-plicit form of this matrix reads Q − = (sin( φ ) σ ph z + cos( φ ) σ ph x ) ⊗ σ RLx ⊗ σ abx . (63)This matrix also satisfies Q − = , and the rela-tions (60,61) stay intact.We close this section with an important statement: theEilenberger functions Q + and Q − of the terminals definethe “Majorana number” of the wire as M = Pf( σ RLz Q + ) Pf( σ RLz Q − ) , (64)in terms of the product of two Z topological invariantsgiven by Eq. (50). It is shown in Appendix C, that for M = − V. NUMERICS
We next turn to the discussion of our numerical resultsobtained for the setup shown in Fig. 1. We have solvedthe quasiclassical equations according to the algorithm ofsection IV E and from there computed the local densityof states (LDoS) ν L ( (cid:15) ) = (2 πv ) − Re tr( σ RLz Q L ) at theleft end of the wire close to the tunnel barrier. In theactual calculations we shifted the energy into the com-plex plane, (cid:15) → (cid:15) + i Γ. This shift accounts for the factthat in the “real” system states may escape to a con-tinuum of lead scattering states, which gives them thestatus of quantum resonances of finite lifetime ∼ Γ − .The corresponding decay rate is given by the standard1golden rule expression, Γ ∼ g T δ where δ ∼ πv/ N L isthe mean level spacing in the scattering region of size L .In principle, one might numerically compute the broad-ening by an extension of the numerical setup so as toinclude and extended normal metallic scattering regionto the left of the tunnel barrier. This, however, wouldslow the performance of the numerics which is why weprefer to introduce the broadening “by hand”. At anyrate, the Majorana peak, present for M = −
1, will ac-quire a finite width ∼ Γ. The DoS profiles obtained as aresult of the procedure outlined above are presented anddiscussed in section II above.
VI. CONCLUSIONS
In this paper we have adapted the Eilenberger quasi-classical approach to the specific conditions of a class Dtopological quantum wire. The most striking feature ofthe quasiclassical approach is that the Green functionsof the system are described by ordinary, and hence eas-ily solvable differential equations. A numerical solutionof these equations for given realizations of the disorderproduces accurate information on the spectral proper-ties of reasonably realistic model systems hosting Ma-jorana boundary states. Our results confirm the predic-tions made in Ref. for a more idealized system, viz. thatdisorder generates spectral peaks which can be confusedfor genuine Majorana peaks. It thus seems that an unam-biguous detection of the Majorana might call for a mea-surement scheme beyond direct tunneling spectroscopy. VII. ACKNOWLEDGMENTS
We thank C. W. J. Beenakker, F. von Oppen andM. Feigel’man for fruitful discussions. This work wassupported by the collaborative research grant SFB/TR12of the Deutsche Forschungsgemeinschaft.
Appendix A: Derivation of Eq. (42)
The goal of this Appendix is to derive the structure(42) by explicit calculation. Our starting point is themomentum representation of the equations of motion forthe quasiclassical Green function g ( x, x (cid:48) ; (cid:15) ) ≡ g ( x − x (cid:48) ; (cid:15) ), (cid:18) − p + iv L (cid:15) (cid:19) g ( p, (cid:15) ) = 1 . (A1)Assuming a diagonal representation as in (39), we obtain g ( p, (cid:15) ) = − T p − i ˆ λv ⊗ σ RLz T − . (A2)The inspection of eigenvalues λ ± shows that for the re-tarded Green’s function, i.e. at (cid:15) → (cid:15) + i
0, one has Re( λ ± ) >
0. This gives g ( x, x (cid:48) ; (cid:15) ) = T (cid:90) dp π e ipx p − i ˆ λv ⊗ σ RLz T − = (A3)= − i T ( σ RLz + sgn( x − x (cid:48) )) e − ˆ λv | x − x (cid:48) | T − . Sending x (cid:48) → x and comparing to (35), we obtain theidentifications (41) and (42). Appendix B: Boundary conditions
In this Appendix we derive the boundary condi-tions (58,59) for the Q -matrix. The method we use isadapted from the circuit theory of Refs. . For sim-plicity, we consider the 4 × g ( x, x (cid:48) ) ≡ g ( x, x (cid:48) ; (cid:15) ) to the right of the right terminal, x, x (cid:48) ≥ x R .By assumption, the generator L (cid:15) (cid:12)(cid:12) x>x R is constant inspace, and this means that the equations (30) defining g admit the solutions g ( x, x R ) = e − L (cid:15) ( x − x R ) g ( x R + 0 , x R ) ,g ( x R , x ) = g ( x R , x R + 0) e L (cid:15) ( x − x R ) , where we have again set v = 1 for notational simplicity.According to Eq. (A3), we have g ( x R + 0 , x R ) = − i Q R + ) ,g ( x R , x R + 0) = − i Q R − ) , (B1)so that g ( x, x R ) = − i e − L (cid:15) ( x − x R ) ( Q R + ) ,g ( x R , x ) = − i Q R + ) e L (cid:15) ( x − x R ) . Defining ¯ g = T − gT and ¯ Q R = T − Q R T , we next trans-form to the representation (39) to arrive at¯ g ( x, x R ) = − i e − ˆ λ ⊗ σ RL ( x − x R ) ( ¯ Q R + ) , ¯ g ( x R , x ) = − i Q R + ) e ˆ λ ⊗ σ RL ( x − x R ) . The key observation now is that the matrix functions¯ g must remain finite as x → ∞ . Due do the posi-tivity of the matrix ˆ λ , this is equivalent to the condi-tion ( − σ RL )( ˆ Q R + ) = ( − σ RL ⊗ )( ˆ Q R + ) =( − Λ)( ˆ Q R + ) = 0. By the same token we have( ¯ Q R + )( + Λ) = 0. Finally, noting that T Λ T − = Q + represents the asymptotic form of the Q -matrix deep inthe superconductor, we arrive at Eqs. (59). Eqs. (58) areshown in an analogous way.2We are now in position to derive Eqs. (60) and (61).To this end, we take Eqs. (58) and multiply it by thetransfer matrix M from the left and by M − from theright. Bearing in mind Eq. (55), one obtains − M Q − M − Q R + M Q − M − − Q + = 0 . (B2)Adding this relation to Eqs. (59), we arrive at2 · − M Q − M − Q R − Q + Q R + M Q − M − − Q + = 0 , (B3)which in one more step yields the result (60). Eq. (61) isproven analogously. Appendix C: Majorana mode
In this Appendix we discuss the “Majorana number” M and show how the Majorana state emerges from thequasiclassical Eilenberger function.We consider the spectrum { E j } of Andreev boundstates which follows from the poles of the Q -function (60,61). If we denote D ( (cid:15) ) = Q + ( (cid:15) ) + M ( (cid:15) ) Q − ( (cid:15) ) M − ( (cid:15) ) , (C1)then the energies E j are solutions of the secular equa-tion det D ( E j ) = 0. The Majorana state, if it exists,corresponds to E = 0.To proceed, we introduce matrices ˜ Q ± ( (cid:15) ) = Q ± ( (cid:15) ) σ RL z ,and the secular matrix ˜ D ( (cid:15) ) = D ( (cid:15) ) σ RL z . Since in the sub-gap interval of energies there is no distinction between theretarded and advanced Green’s function, the unitarity re-lation ( G R/A ) † = G A/R , the basic definitions Eqs. (28)and (35) imply ˜ Q †± ( (cid:15) ) = − ˜ Q ± ( (cid:15) ). Further, particle-holesymmetry (37) yields ˜ Q T ± ( − (cid:15) ) = − ˜ Q ± ( (cid:15) ). Taking intoaccount the class D symmetry of the transfer matrix,Eq. (57), the secular matrix takes a form˜ D ( (cid:15) ) = ˜ Q + ( (cid:15) ) + M ( (cid:15) ) ˜ Q − ( (cid:15) ) M T ( − (cid:15) ) . (C2)It satisfies ˜ D T ( (cid:15) ) = − ˜ D ( − (cid:15) ), and thereby guarantees thatAndreev bound states appear in pairs ± E j .One may ask now whether a zero energy solution ofthe secular equation exists or not. At (cid:15) = 0, the particlehole symmetry puts tighter restrictions on the matrices˜ Q ± , M and ˜ D (we omit the energy argument for brevity).One gets ˜ Q T ± = − ˜ Q ± , ˜ Q ∗± = ˜ Q ± , (C3) σ RL z M T σ RL z = M − , M ∗ = M, (C4)˜ D T = − ˜ D , ˜ D ∗ = ˜ D . (C5)We thus observe that both ˜ Q ± and ˜ D are real antisym-metric matrices of size 8 N × N , which enables us torewrite the secular equation in terms of a PfaffianDet ˜ D = (cid:104) Pf( ˜ D ) (cid:105) = (cid:104) Pf( ˜ Q + + M ˜ Q − M T ) (cid:105) = 0 . (C6) Let us denote by Ω L the left null space of the matrix ˜ D ,i.e. any bra (cid:104) φ | ∈ Ω L by definition satisfies (cid:104) φ | ˜ D = 0.Since ˜ D is antisymmetric, the dimension of its null spaceis even, dim Ω L = 2 N . At this stage we make use of amathematical lemma proven in Appendix A of Ref. 24:the parity of the number N can be expressed as( − N = Pf ˜ Q + Pf [ M ˜ Q − M T ] = Det M Pf ˜ Q − Pf ˜ Q + . (C7)The particle-hole symmetry (C4) implies that Det M = ± M ( x, x (cid:48) ) isa continuous function of its arguments, with initial value M ( x = x (cid:48) ) = . We thus conclude that Det M = +1,so that the parity of N is determined by the terminalconfigurations, ( − N = Pf ˜ Q − Pf ˜ Q + . (C8)This parity is equal to the “Majorana number” M intro-duced in section IV E.Let us now focus on the most interesting case N =1, corresponding, as we will see, to a single Majoranamode . In this case we have dim Ω L = 2, and the nullspace is spanned by two linearly independent vectors. Let (cid:104) φ | ∈ Ω L be the first basis vector. It is easy to check that (cid:104) φ | = (cid:104) φ | Q + ∈ Ω L , and we may choose this state forthe second basis vector in Ω L . Indeed, for any (cid:104) φ | ∈ Ω L one has (cid:104) φ | ˜ D = (cid:104) φ |D = 0. Using the definition of D ,Eq. (C1), we deduce that (cid:104) φ | Q + M = −(cid:104) φ | M Q − . (C9)This relation enables us to evaluate (cid:104) φ |D as (cid:104) φ |D = (cid:104) φ | (cid:0) + Q + M Q − M − (cid:1) = (cid:104) φ | + ( (cid:104) φ | Q + M ) Q − M − = (cid:104) φ | − (cid:104) φ | M Q − Q − M − = 0 , (C10)and hence we proved that (cid:104) φ | ∈ Ω L . Since D is areal antisymmetric matrix, its right null space is ob-tained as Ω R = (Ω L ) † . In other words, the two kets | φ , (cid:105) = ( (cid:104) φ , ) | ) † satisfy the relation ˜ D| φ , (cid:105) = 0.Let us look at the Eilenberger function Q R ( (cid:15) ) aroundits pole at (cid:15) = 0. According to Eq. (60), it can be repre-sented by two equivalent equations, Q R = + 2 D − ( − Q + ) , (C11) Q R = − + 2( + Q + ) D − . (C12)Multiplication by σ RL3 yields σ RL3 ˜ Q R σ RL3 = σ RL3 + 2 ˜ D − ( − Q + ) , (C13) σ RL3 ˜ Q R σ RL3 = − σ RL3 + 2( − Q T+ ) ˜ D − . (C14)At (cid:15) → D − ∼ R /(cid:15) , where the matrix R is its residue at zero energy.At this stage it is advantageous to introduce a new bra (cid:104) χ ± | = (cid:104) φ | ± (cid:104) φ | = (cid:104) φ | ( ± Q + ) , (C15)3and ket | χ ± (cid:105) = ( ± Q T+ ) | φ (cid:105) , (C16)basis in the null spaces Ω L and Ω R , resp. The bra basis (cid:104) χ ± | is not orthogonal and we can denote by | η ± (cid:105) the ketbasis, which is dual to it, (cid:104) χ σ | η σ (cid:48) (cid:105) = δ σσ (cid:48) . Using thesedefinitions, we may formulate resolutions of unity, = | η σ (cid:105)(cid:104) χ σ | = | χ σ (cid:105)(cid:104) η σ | , where σ = ± is summed over.Now let us look at the singular part of the matrix ˜ Q R , σ RLz ˜ Q sing R σ RLz = 2 (cid:15) R ( − Q + ) , = 2 (cid:15) ( − Q T+ ) R . (C17)We note that the matrix P − = ( − Q + ) acts as theprojector in the space Ω L , since (cid:104) χ + | P − = 0 , and (cid:104) χ − | P − = (cid:104) χ − | . (C18)Similar relations hold in the ket space Ω R : P T − | χ + (cid:105) = 0 , and P T − | χ − (cid:105) = | χ − (cid:105) . (C19)Equivalently, we may write P − = | η − (cid:105)(cid:104) χ − | , P T − = | χ − (cid:105)(cid:104) η − | . Using these properties, it follows σ RLz ˜ Q sing R σ RLz = (cid:15) R P − = (cid:15) P T − R . Multiplying these relations by P − fromthe right, and using the projector property P − = P − , weobtain the relation σ RLz ˜ Q sing R σ RLz = 4 (cid:15) P T − R P − = 4 (cid:15) | χ − (cid:105)(cid:104) η − |R| η − (cid:105)(cid:104) χ − | , or ˜ Q sing R = 4 R −− (cid:15) σ RLz | χ − (cid:105)(cid:104) χ − | σ RLz , where we defined R −− ≡ (cid:104) η − |R| η − (cid:105) . Taking into ac-count our definition of the quasiclassical Green’s func-tion, Eqs. (28) and (35), we finally obtain that thesingular part of the propagator around zero energy at x = x (cid:48) = x R takes the form G ( x R , x R ; (cid:15) ) ∼ (cid:18) R −− i(cid:15) (cid:19) ˆ v − / σ RLz | χ − (cid:105)(cid:104) χ − | σ RLz ˆ v − / , (C20)where ˆ v is the diagonal velocity matrix.This expression can be compared with the spectral de-composition of the Green’s function, G ( x, x (cid:48) , (cid:15) ) = | ψ ( x ) (cid:105)(cid:104) ψ ( x (cid:48) ) | (cid:15) + (cid:88) j | ψ E j ( x ) (cid:105)(cid:104) ψ E j ( x (cid:48) ) | (cid:15) − E j + (cid:32)(cid:90) + ∞ E g d E + (cid:90) − E g −∞ d E (cid:33) | ψ E ( x ) (cid:105)(cid:104) ψ E ( x (cid:48) ) | (cid:15) − E , where | ψ E ( x ) (cid:105) are normalized eigenfunctions of the BdGHamiltonian, E g is the gap in the spectrum of the wire,the sum is going over the set of subgap Andreev levels(we have particle-hole symmetry E − j = − E j ) and thefirst term is present if the system contains a Majoranastate. One thus concludes that (cid:18) R −− i ˆ v (cid:19) − / σ RLz | χ − (cid:105) = | ψ ( x R ) (cid:105) (C21)is the amplitude of the Majorana particle at the point x R . The amplitude of the Majorana state at any otherpoint x can be obtained from | ψ ( x R ) (cid:105) by applying thetransfer matrix M ( x, x R ).To conclude, we have shown that if the Majorananumber of our system is topologically non-trivial, i.e. M = −
1, then the Green’s function has the pole at zeroenergy. Up to a prefactor, the residue at zero energy isthen the projector on the one dimensional linear subspacespanned by the Majorana particle.
Appendix D: Matrix Q − of a tunnel junction In this appendix we show that the boundary condi-tions for the Eilenberger matrix Q of a spin-orbit wireterminated at the left end by a (infinitely high) tunnelbarrier are equivalent to algebraic relations (58) with theeffective matrix Q − given by Eq. (63).Let us work in the original particle-hole basis wherethe spinor Ψ takes the form specified by Eq. (14). Westart by rewriting the left boundary conditions (62) in amatrix form. If one defines the reflection matrixˆ r = (cid:18) e iφ σ RLx e − iφ σ RLx (cid:19) (D1)in the particle space, then ˆ r ∗ is the reflection matrix forholes. Consequently, the spinor Ψ( x L ) satisfies the con-dition (cid:16) − ˆ R (cid:17) Ψ( x L ) = 0 , ˆ R = (cid:18) ˆ r
00 ˆ r ∗ (cid:19) . (D2)The full reflection matrix ˆ R is unitary obeying theparticle-hole symmetry, σ ph x R T σ ph x = R . Hence theboundary condition for the bar spinor ¯Ψ = ( σ ph x Ψ) T takesthe similar form ¯Ψ( x L ) (cid:16) − ˆ R (cid:17) = 0 . (D3)According to the definition (28), the Green’s function g e ( x, x (cid:48) ) inherits the boundary condition right to thetunnel junction from those of the direct product of twospinors, Ψ( x ) ⊗ (cid:0) ¯Ψ( x (cid:48) ) σ RLz (cid:1) . We thus conclude that (cid:0) − ˆ R (cid:1) g ( x L , x ; (cid:15) ) = 0 ,g ( x, x L ; (cid:15) ) (cid:0) + ˆ R (cid:1) = 0 . σ RLz ˆ r σ RLz = − ˆ r . The above conditions arevalid for any x > x L . Sending now x → x L + 0 and usingthe relations (B1) written for the matrix Q L , one obtains( − ˆ R )( − Q L ) = 0 , ( + Q L )( − ˆ R ) = 0 . We see that these boundary conditions are equivalent tothe algebraic conditions (58) if one identifies Q − = − ˆ R . The normalization ˆ R = guaranties that Q − belongs tothe manifold of the quasiclassical Eilenberger functions.Transforming matrix Q − into the Majorana representa-tion (20) we finally obtain the result (63).The scattering phase φ is the parameter of the matrix Q − which may depend on the band index n and charac-terizes the tunnel junction. In our numerical simulationswe have used φ = π for all bands, which corresponds tothe infinitely high barrier as compared to the energies ofAndreev bound states. R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev.Lett. , 077001 (2010). Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. , 177002 (2010). V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P.A. M. Bakkers, and L. P. Kouwenhoven, Science ,1003 (2012). M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff,and H. Q. Xu, arXiv:1204.4130 (2012). A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, andH. Shtrikman, arXiv:1205.7073 (2012). D. Bagrets and A. Altland, Phys. Rev. Lett. , 227005(2012). J. Liu, A. C. Potter, K. T. Law, and P. A. Lee, Phys. Rev.Lett. , 267002 (2012). C. W. J. Beenakker, D. I. Pikulin, T. Hyart, and J. P.Dahlhaus, New J. Phys. , 125011 (2012). G. Eilenbeger, Zeitschrift fur Physik , 195 (1968). A. Larkin and Y. Ovchinnikov, Sov. Phys. JETP , 1200(1969). J. Rammer and H. Smith,
Quantum Field-theoretical Meth-ods in Transport Theory of Metals , Reviews of modernphysics (American Physical Society, 1986). M. Eschrig, Phys. Rev. B , 9061 (2000). A. Altland, B. D. Simons, and D. Taras-Semchuk, Ad-vances in Physics , 321 (2000). M. V. Feigel’man, A. I. Larkin, and M. A. Skvortsov,Phys. Rev. B , 12361 (2000). Y. Nazarov and Y. Blanter,
Quantum Transport: Introduc-tion to Nanoscience (Cambridge University Press, 2009). A. Altland and M. R. Zirnbauer, Phys. Rev. B , 1142(1997). A refined variant of Eq. (1) has recently been derived inRef. 26. It was found that for extremely low temperatures, T (cid:28) Γ the tunneling current may contain a dip reflectingthe mutual cancellation of electron and hole current con-tributions. Our discussion is thus tacitly assumes T (cid:38) Γ. R. M. Lutchyn, T. D. Stanescu, and S. Das Sarma, Phys.Rev. Lett. , 127001 (2011). A. Shelankov and M. Ozana, Phys. Rev. B , 7077 (2000). Y. V. Nazarov, Superlattices and Microstructures , 1221(1999). A. Kitaev, Phys. Usp. , 131 (2001). M. Wimmer, A. R. Akhmerov, J. P. Dahlhaus, andC. W. J. Beenakker, New J. Phys. , 053016 (2011). A. Cottet, D. Huertas-Hernando, W. Belzig, and Y. V.Nazarov, Phys. Rev. B , 184511 (2009). I. C. Fulga, F. Hassler, A. R. Akhmerov, and C. W. J.Beenakker, Phys. Rev. B , 155429 (2011). In case of odd
N ≥
N ≥ N -fold degeneracy of a zero level is accidental and nottopologically protected. It can be reduced down to N = 1or 0 by continuous distortion of the transfer matrix M .26