Quasiclassical theory of disordered multi-channel Majorana quantum wires

Multi-channel spin-orbit quantum wires, when subjected to a magnetic field and proximity coupled to s-wave superconductor, may support Majorana states. We study what happens to these systems in the presence of disorder. Inspired by the widely established theoretical methods of mesoscopic superconductivity, we develop a la Eilenberger a quasiclassical approach to topological nanowires valid in the limit of strong spin-orbit coupling. We find that the"Majorana number", distinguishing between the state with Majorana fermion (symmetry class B) and no Majorana (class D), is given by the product of two Pfaffians of gapped quasiclassical Green's functions fixed by right and left terminals connected to the wire. A numerical solution of the Eilenberger equations reveals that 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 the Majorana peak. In this way we confirm the robustness of our previous conclusions [Phys. Rev. Lett. 109, 227005 (2012)] on a more restrictive system setup. Generally speaking, we find that the quasiclassical approach provides a highly efficient means to address disordered class D superconductors both in the presence and absence of topological structures.


I. INTRODUCTION
Semiconductor quantum wires proximity coupled to a conventional superconductor and subject to a magnetic field may support Majorana fermion edge states 1,2 . Building on relatively conventional device technology, this proposed realization of the otherwise evasive Majorana fermion has triggered a wave of theoretical and experimental activity, which culminated in the recent report of a successful observation by several experimental groups [3][4][5] . In these experiments, evidence for the presence of a Majorana is drawn from the observation of a zero bias peak in the tunneling conductance into the wire. While the observed signal appears to be naturally explained in terms of a Majorana end state, two of us have pointed out that a midgap peak might be generated by an unrelated mechanism 6 : in the presence of even very moderate amounts of disorder, the semiconductor wire supports a zero energy "spectral peak" (an accumulation of spectral weight at zero energy) which resembles the Majorana peak in practically all relevant aspects. Specifically, it is (i) rigidly locked to zero energy, (ii) is of narrow width of O(δ), where δ is the single particle level spacing, (iii) carries integrated spectral weight O(1), and (iv) relies on parametric conditions (with regard to spin orbit interaction, proximity coupling, magnetic field, and chemical potential) identical to those required by Majorana state formation. What makes the spectral peak distinct from the Majorana is that it relies on the presence of a moderate amount of disorder, viz. impurity scattering strong enough to couple neighboring single particle Andreev levels. Besides, the spectral peak is vulnerable to temperature induced dephasing. While this marks a difference to the robust Majorana state, the reported experimental data does show strong sensitivity to temperature, which may either be due to an intrinsic sensitivity of the peak, or due to a temperature induced diminishing of the measurement sensitivity, or both. In either case, the situation looks inconclusive in this regard.
Generally speaking, the results of Ref. 6, as well as those of Refs. 7 and 8 suggest that the observation of a midgap anomaly in the tunneling conductance might be due to either mechanism, disorder peak, Majorana peak, or a superposition of the two, and this calls for further research.
Our previous study was based on an analytically tractable idealization of a semiconductor quantum wire subject to a magnetic field sweep. In the present paper, we will explore the role of disorder scattering within a model much closer to the experimental setup. The price to be payed for this more realistic description is that a fully analytic treatment is out of the question. Instead, we will employ a semi-analytic approach based on the formalism of quasiclassical Green functions. Introduced in the late sixties 9,10 , the latter has become an indispensable tool in mesoscopic superconductivity [11][12][13][14] and quantum transport in general 15 . We here argue that quasiclassical methods are, in fact, tailor made to the modeling of Majorana quantum wires (or, more generally, quasi onedimensional topological superconductors.) Specifically, we will show that in the problem at hand, the quasiclassical "approximation" is actually very mild. Further, the quasiclassical Green's function affords a convenient description of the topological signatures of the system in terms of Pfaffians. Finally, the approach can be applied to systems for a given realization of the disorder, and at numerical cost much lower than that of exact diagonalization approaches. As a result, we will be in a position to accurately describe local spectral properties within a reasonably realistic model of a topological multi-channel superconductor. As we are going to discuss below, our findings support the principal statement made in Ref. 6.
The rest of the paper is organized as follows. In section II, we discuss the principal role played by disorder in the system, the idea of the quasiclassical approach, and its main results. In section III we specify our model system, the quasi-classical approach is introduced in section IV, and in V we discuss the numerical solution of the quasiclassical equations. We conclude in section VI. A number of technical details are relegated to Appendices.

II. QUALTIATIVE DISCUSSION AND RESULTS
A schematic of the device currently under experimental investigation is shown in Fig. 1. A semiconductor quantum wire subjected to strong spin orbit interaction is brought in contact to a superconductor (S), and, via a tunnel barrier (T) to a normal metal lead (N). The application of a small excess voltage, V , to the latter induces a tunnel current into the central region. The differential conductance dI/dV probes the (tunneling) density of states at an energy V (units e = = c = 1 throughout) relative to the systems chemical potential, µ. The physics 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, one contribution to the band center density of states is provided by a Majorana bound state localized at the tunnel barrier. The second, spectral contribution is generated by a conspiracy all other low lying quasiparticle states in the system. Technically, these are Andreev bound states forming at energies ±E j in the region between the tunnel barrier and the superconductor. The number of these states increases with the extension of the wire -a few hundred nanometers in the experiment -and the number of transverse channels below the chemical potential. In the presence of disorder, the entity of these states defines an effective "quantum dot". The proximity of a superconductor, and the breaking of both spin rotation and time reversal symmetry imply that the system belongs to symmetry class D 16 .
The symmetry of class D random systems implies a clustering of levels at zero energy. Loosely speaking, the conventional level repulsion of random spectra turns into a zero energy level attraction. On a resolution limited to scales of order of the mean level spacing, the zero energy density of states (DoS) is enhanced by a factor of two relative to the mean background, i.e. it shows a peak. This phenomenon manifests itself at relatively small sample-to-sample fluctuations, i.e. the peak is a sample specific effect. In passing we note that the weak anti-localization phenomenon discussed in Ref. 8 rests on the same principal mechanism of midgap quantum interference.
Below, we will explore the phenomenon of spectral peak formation in a setting modelled to closely mimic the "experimental reality". To be more specific, we consider a semiconductor wire supporting a number of N > 1 transverse channels below the chemical potential. We assume a value of the chemical potential such that the highest lying of these channels is "topological" (chemical potential falling into the gap opened by the simultaneous presence of order parameter and magnetic field.) To quantitatively describe this phenomenon, we generalize the quasiclassical Green function approach to the present setting. Indeed, the phenomena we are interested in manifest themselves on length scales large in comparison to the Fermi wavelength, yet smaller than the relevant coherence length, and this makes them suitable to quasiclassical treatment. The quasiclassical theory of the present paper is formulated in terms of the Eilenberger functionQ(x; ), which is a position and energy dependent matrix of size 8N × 8N . The factor of 8 = 2 × 2 × 2 accounts for spin, chiral (left/right modes) and particlehole degrees of freedom. It will turn out that the structure of the theory is most transparently exposed in the so-called Majorana basis, where the Eilenberger function becomes a real antisymmetric matrix. The Pfaffian of that matrix, Pf(Q) will be seen to assume two values (±1), which locally (in space) identify the Z 2 invariant of the underlying one-dimensional class D superconductor, in dependence on system parameters: for values of the magnetic field smaller or larger than a critical field, B c ≡ ∆ 2 + µ 2 , where µ and ∆ are chemical potential and bulk order parameter, the invariant is trivial or non-trivial, resp. In the latter case, the wire supports a Majorana state at the barrier, in the former it doesn't. Fig. 2 shows a profile of the local DoS (LDoS) computed at the left end of the wire for typical system parameters as detailed below. The green dashed and red solid curve correspond to a situation without and with Majorana state. Within our model, the states in the wire carry a narrow width, ∼ g T δ, reflecting the possibility of decay through the tunnel barrier into the adjacent lead. (Here, δ is the mean level spacing of the wire, and g T 1 the barrier tunneling conductance.) This broadening accounts for the finite width of the Majorana peak in the topological regime. Loosely speaking, the negative shoulders observable next to the center peak are due to the repulsion of adjacent levels of the center Majorana level. A more substantial explanation is as follows: for an odd parity of the total number of levels -a signature of the topological regimes -the disordered quantum system falls into symmetry class B, rather than D. (Class B is the designation for a system with the same symmetries as a class D system, yet odd number of levels.) A universal signature of class B is a negative spectral peak at zero energy (the negative of the positive class D peak), superimposed on a single δ-peak (the Majorana). The joint signature of these two structures is seen in the solid curve in Fig. 2.
At resolutions limited to values ∼ δ, the superposition of the Majorana and the class B peak looks next to indistinguishable from the class D peak (dashed), and this similarity of unrelated structures might interfere with the unambiguous observation of the Majorana by tunneling spectroscopy. Indeed, the differential tunneling conductance monitored in experiment,  The profiles of the curves shown in Fig. 2 were computed for a two-channel quantum wire at a mean free path l L of the order of the system size. We are addressing a system at the interface between the ballistic and the localized regime. In view of these system parameters it is remarkable that the DoS profiles in Fig. 2 show striking similarity to the average DoS of a class D and B random matrix model 16 . For comparison the average DoS of a class B and D random matrix Hamiltonian is shown as an inset in Fig. 2. The similarity of the results indicates that the system of subgap states in our system behaves as if it formed an effective chaotic quantum dot localized in the vicinity of the left system boundary (right to the tunnel barrier). Taking into account spin, chiral and channel quantum numbers, the mean level spacing in such dot is given by δ πv/2N L. In the presence of magnetic field the BCS gap in the superconducting region of the wire reads − = |B − ∆ 2 + µ 2 |. The number of subgap Andreev levels forming the effective dot is thus given by N levels 2 − /δ.
The profiles shown in Fig average ν 2 L . The relatively minor deviation between average and typical DoS demonstrates that the standard deviation characterizing the strength of mesoscopic fluctuations is relatively small.
The weakness of fluctuations implies that the disorder peak is a realization specific phenomenon. This is demonstrated explicitly in Fig. 4, where two unaveraged DoS profiles individual disorder configurations are shown. The data is for a "non-topological" system, no non-degenerate zero energy level is present. Depending on whether the pair of lowest lying Andreev states (+ min , − min ) exceeds the level broadening Γ, the DoS enhancement will assume the form of a single peak (magenta solid), or a split peak (green dashed). In either case, an excess DoS of integrated spectral weight ∼ 1 is present.
We finally note that both the Majorana peak and the D spectral peak crucially rely on the presence of a magnetic field. In the absence of a field time-reversal symmetry is restored, and the wire turns into a member of symmetry class DIII. Such systems display a DoS depletion at zero energy, rather than a peak. In other words, the spectral peak discussed here will disappear along with the Majo-rana resonance.
In the rest of the paper we discuss how the results summarized here were obtained by a combination of quasiclassical and numerical methods.

III. MODEL OF THE MAJORANA NANOWIRE
We consider a multi-band quantum wire of width L z , subject to Rashba spin-orbit coupling, proximity coupling to an s-wave superconductor, and to a magnetic field 18 . We choose coordinates such that the wire lies along the x-axis, parallel to the magnetic field, the y-axis is perpendicular to the surface of the superconductor, and the spin orbit field is pointing along the z-axis.
In the rest of this section we will specify the Bogoliubov-de-Gennes (BdG) Hamiltonian describing this system in the presence of disorder (section III A). We will then linearize the electron spectrum in the limit of strong spin-orbit coupling thereby introducing a description in terms of right and left one-dimensional chiral fermions (section III B), and finally transform the Hamiltonian into Majorana basis (section III C). Experts in the description of topological quantum wires may proceed directly to the end of the section.

A. Bogoliubov-de-Gennes Hamiltonian
Introducing a four component spinor Ψ = (ψ ↑ , ψ ↓ ,ψ ↑ ,ψ ↓ ) in the product of spin and particlehole spaces, the BdG Hamiltonian H describing the system in the xz-plane reads where H is the external field, we have taken into account the transverse momentum (−i∂ z ) in the Rashba term, andŴ (x) is the random disorder Hamiltonian. The Pauli matricesŝ operate in spin space.
We proceed by introducing a system of transverse wave functions, {Φ σ n (z)}, which leads to a linear decomposition of the Grassmann fields (σ =↑, ↓) as Defining the Hamiltonian then takes the form where the one-dimensional Hamitonianĥ In the case of an ideal waveguide the transverse wavefunctions are Φ σ n (z) = 2/L z sin(nπz/L z ), so that µ (n) z = π 2 n 2 /(2mL 2 z ) and the matrix elements of the spin-orbit interaction read (8) Let us now assume a thin wire, L z l so = /(mα), where l so is the spin-orbit length. In this case the ma- can be treated as perturbations. To this end we introduce a unitary transforma-tionÛ, which brings the high energy part of the Hamiltonian, (µ nm ), to diagonal form. Due to the weakness of the perturbation, the transformation is close to unity,Û = exp(iX) 1 + iX, with generatorŝ X nm = O(L z /l so ). To first-order perturbation theory their explicit form readŝ The transformation Ψ →ÛΨ andΨ →ΨÛ † generates the approximate Hamiltonian Here we have redefined our notation for the disorder potentialŴ →ÛŴÛ † , and introduced a small correction to the quasiparticle Hamiltonian,V = α [X,ŝ z ] ∂ x , and to the order parameter, δ∆ = −∆[X,ŝ y ]. One can estimate the matrix elements of these operators as V nm ∼ (L z /l so ) and (δ∆) nm ∼ ∆(L z /l so ), where is the quasiparticle energy. Since we are primary interested in the effects of disorder, we will limit our discussion to the situation when the low-energy physics is disorder dominated, i.e. the random potentialŴ > max{V , δ∆} masks the off-diagonal matrix elements of the deterministic Hamiltonian. Having this in mind we thus omit botĥ V and δ∆ throughout the paper. However, this assumption does not affect our results in qualitative ways.
Assume that the chemical potential lies close to the bottom of the N -th band, µ µ (N ) z . We refer to this band as the "topological band", since in the absence of interband scattering it defines whether the quantum wire is in the topological phase or not. The condition of the topologically non-trivial phase reads We also assume a hierarchy of energy scales µ where E so = mα 2 /2 is the spin-orbit energy. This condition implies that all other channels with the band index n < N are in the trivial phase, since for them In the current experiments E so ∆, i.e. our theory is on the border of applicability (see, however, the discussion in section IV).

B. One-dimensional chiral fermions
For strong spin-orbit coupling, E so B, each band is characterized by two Fermi momenta, as shown in Fig. 5. We aim to construct a low energy Hamiltonian describing the system at energy scales E so . To this end we represent the fields ψ around the Fermi momenta. Modes with Fermi momenta of equal modulus define a conduction channel. We observe that each channel belongs to one of the two subsets, a or b. In the a-channel R-movers have spin up, and Lmovers have spin down, while for the channels of type b spins are reversed. The a-channel of band index N plays a distinguished role in that it has zero Fermi momentum, k a N = 0. This particular mode defines what we call the "topological channel". For max{B, ∆} E so it is the only 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 be safely 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 implicit. The spinor Ψ acts in an 8 × N -dimensional space, defined by the direct product of band (index n), channel (a/b), chiral (R/L) and particle-hole spaces. In this representation the low energy Hamiltonian is given by whereĥ Here is the Fermi velocity of the n-th channel, the chemical potential µ is now defined relative to the energy µ (N ) z , andŴ nm (x) is a random 8 × 8 matrix satisfying the hermiticity condition (Ŵ nm ) † =Ŵ mn . In deriving this Hamiltonian we have neglected oscillatory terms with phases ix(±k a ±k b ). This is justified because the typical lengths involved in the problem (ξ ∼ ∆/v and l B ∼ B/α) exceed by far the Fermi length λ F ∼ max(1/k a , 1/k b ) determined by the large spin-orbit energy E so . We have also used that the order parameter ∆(x) can be chosen to be real.

C. Majorana representation
The 8N × 8N first-quantized matrix defining the Hamiltonian (15) satisfies to the p-h symmetry, which is the defining condition for a class D Hamiltonian (here, the transposition acts on kinetic term as However, all what follows will be more con-veniently formulated in an alternative representation, in which the symmetry assumes a different form: for each band n we define a set of eight Majorana fields (and analogous relations for the L-movers, with spin reversed), which we combine into the 8N -spinor The two spinors χ and Ψ are related by the unitary transformationχ We combine the Majorana fields of a and b channels as ξ R = (ξ R a , ξ R b ) T , and reorder the spinor components as, In this representation (which will be used throughout the rest of the paper) the Hamiltonian (15) takes the form where the deterministic part readŝ and the random matrices are constructed aŝ in terms of real (symmetric) and imaginary (antisymmetric) parts of the random matrices (Ŵ ) RL =ŵ RL 1 + iŵ RL 2 etc. The remaining blocks of theŴ -matrix are defined byŴ In the Majorana basis (20), the class D symmetry is expressed through the antisymmetryĤ T = −Ĥ.
We finally specify the statistics of disorder. We choosê W (x) to be a δ-correlated and Gaussian distributed random matrix of size 8N × 8N with a zero mean value Ŵ (x) and variance Here, the composite indices (i, j etc.) label states in the direct product of band, channel and chiral spaces. The scattering matrices defined in this way break timereversal and spin rotation symmetry (e.g., random spinflip scattering caused by random spin-orbit terms is included in (26).) The strength of disorder set by the coefficient γ w translates into the "golden rule" scattering rate of the normal conducting (i.e. superconductor decoupled) quantum wire.

IV. QUASICLASSICAL APPROACH
The kinetic term in the low energy Hamiltonian (21) which was discussed in the previous section is linear in momentum, and this facilitates the formulation of quasiclassical equations of motion (aka Eilenberger equations) for the model at hand 9 . We here review the construction of these equations in a manner closely following the spirit of Refs. 19 and 20. After introducing the basics of the method (section IV A) we construct the Eilenberger Q-function in the limit of a single clean "topological channel" (section IV B) and discuss the resulting density of states (section IV C). In section IV D we define the Z 2 topological invariant in terms of the Q-matrix. Section IV E outlines the general construction of the solution to the Eilenberger equation in the inhomogeneous disordered wire with boundary conditions.

A. Eilenberger method
We start by defining, where G R/A (x, x ) ≡ x|( ±i0−H) −1 |x are the retarded and advanced Green's functions of the system. Under transposition (which in our current representation represents the particle-hole symmetry) the function g behaves as i.e. advanced and retarded Green's functions get interchanged. It is not hard to derive the two mutually adjoint differential equations describing the dynamical evolution of g. Here, where matrixω has elements and the operator P is related to the Hamiltonian matrix H as Due to the antisymmetry ofH = −H T , the operator P obeys the particle-hole symmetry We next define the Eilenberger function as where the subtraction of the sgn-function regularizes a discontinuity arising in g at x = x due to the combination of linear derivatives and δ-function inhomogeneity in Eqs. (30). Subtracting the two equations in (30), we then obtain the Eilenberger equation of motion The Eilenberger function Q obeys the particle-hole symmetry and the normalization condition Q 2 (x) = 1, where 1 is the unit matrix (the latter condition can be checked by verifying that Eq. (36) preserve the normalization Q 2 = const.). The unit-value of the normalization constant is fixed by the jump height of the sgn-function in (35). We finally note that the operator P can be straightforwardly constructed 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 approach to the limit (Ŵ = 0) of an infinite clean quantum wire subject to constant B, ∆. In this system all channels are decoupled, and we may concentrate on the 4 × 4 matrix Q (B, µ) describing the "topological" channel a with n = N . (The Eilenberger function of the other channels may be obtained by setting B = 0, rescaling the velocity and transforming ∆ → −∆ in the case of b-type channels. ) We start by introducing the 4 × 4 operator as the reduction of the general operator L to a single channel. The solution Q is then determined by the relations [Q , L ] = 0 and Q 2 = 1. To solve these equations, we assume L to be diagonalized as where the exact form of the (non-unitary) transformation matrix T will not be needed and are the eigenvalues. The defining equations for Q are then solved by matrices of the form where Λ = (±1, . . . , ±1) is a diagonal 4 × 4 matrix containing unit-modular entries in arbitrary configuration. The proper sign structure is determined by causality, i.e. the sign of the infinitesimal offset → ± i0 in the retarded/advanced Green function. That increment enters in the combination ( ±i0)σ RL z , which means that the appropriate matrix structure of the retarded Green's function (opposite for advanced) is given by A more explicit derivation of this structure is detailed in Appendix A.

C. Clean density of states
The density of states in the bulk of the topological wire is given by ν( ) = (2πv) −1 Re tr σ RL z Q . The matrices T diagonalizing Q do not commute with σ RL z , which means that a little extra work is required to evaluate the trace. We start from the representation where P + and P − are projectors on the space of Leigenstates with eigenvalues ±λ + and ±λ − , resp.: That this representation faithfully represents the matrix Q is checked by application of (43) in the eigenbasis where all matrices assume a diagonal form. It remains to obtain a representation of P ± which does not make explicit reference to the diagonalizing matrices T . To this end, notice that (eigenrepresentation understood) L 2 = diag(λ 2 + 1 2 , λ 2 − 1 2 ) = λ 0 1 4 + λP + − λP − . This equation can straightforwardly solved as Substituting this expression into (43), we obtain a representation of Q which makes reference only to the operator (38), and its eigenvalues. Computing the trace, we obtain DoS profiles as shown in Fig. 6. Before discussing the structure of these results, a general remark may be in order: the DoS of a one dimensional quantum system is determined 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, and this matrix seems to be oblivious to the kinetic energy. A closer look, however, shows that information on the band dispersion sneaks in via the nonlinear constraint Q 2 = 1. Indeed, the conservation of the constraint, and the unit value of the normalization are consequences of the linearity of the derivative operator in (30), which in this way co-determines the structure of Q. Inspection of Eq. (43) shows that the DoS contains singularities at the zeros λ ± ( ) = 0, which are located at Let us assume that B > ∆. Fig. 6(a) shows the ensuing DoS profile, along with the underlying dispersion relation for a value of the chemical potential µ < µ c , where defines a critical value where the lower of the DoS singularities, − , 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 the lower band − (k) develops an extremum at finite values of k which manifests in a third van-Hove singularity at the energy 6. Sketch of the density of states and the corresponding dispersion relations in the clean nanowire. The gap is closed at the transition point µ = µc (a), and opened again (b,c). For µ > µ * two minima develop (d).

D. Z2 topological invariant
The symmetry Eq. (37) implies that at zero energy the product σ RL z Q =0 is an antisymmetric 4 × 4 matrix, which implies the existence of a Pfaffian. Due to the signature of Λ the determinant of Q is unity, the same is true for the determinant of the 4 × 4 matrix σ RL z = σ RL z ⊗ 1 2 . Consequently, the Pfaffian of σ RL z Q =0which squares to the determinant of that matrix -may take one of two values, ±1. This motivates the definition of the topological index 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 =0 becomes singular and the index cannot be defined.)

E. Eilenberger equation with disorder
In this section we discuss the formal solution of the Eilenberger equation in the presence of disorder. The solution is "formal" in the sense that the Eilenberger Green function will be a functional of a given realization of the disorder. To obtain practically useful information, one will want to average over different realizations, and this step of the computation needs to be done numerically, as discussed in the next section.
To start with, consider the prototypical system geometry shown in Fig. 7. The terminals indicated at the left and right represent superconducting regions, assumed non-disordered for simplicity. (This is an inconsequential assumption provided the rate of disorder scattering τ −1 ∼ N γ 2 w 1/v n does not exceed the energy gaps ( − or 0 ) in the terminals.) In these regions, the Eilenberger equation can be solved analytically, as discussed in the previous section. Describing the disorder present in the center region in terms of a generalized variant of the quasiclassical evolution operator L, we will show how the left and right asymptotic configuration of the Green function get connected by a transfer matrix, M , functionally depending on the disorder configuration. The ensuing generalized Green function will then be the starting point for our numerical analysis.
To be more specific, we consider a quantum wire where the gap ∆(x) and/or chemical potential µ(x) vary in space in the region |x| < L/2 and saturate to some constants ∆ L/R and µ R/L at x −L/2 or x L/2, respectively (Fig. 7). These constants set asymptotic values of the Q-matrix, where Q ± are constructed using the results of section IV B for the homogeneous profile of ∆, B and µ.
The boundary Green's function Q − and Q + may describe different or equivalent topological phases of the wire. We denote the Q-matrices obtained at the interface between the asymptotic superconducting regions, and the center disordered region, resp., as where x R = −x L = L/2. These two configurations are related by a transfer matrix, For arbitrary positions x and x the formal expression for the transfer matrix M (x, x ) at given energy follows 7. Disordered spin-orbit wire connected to two ideal superconducting terminals which are described by the Eilenberger functions Q+ and Q−. The Q-matrix at the boundaries between the scattering region and terminals is denoted by QR and QL. In the superconductors Q-matrix rapidly converges to either Q+ or Q− on a scale of coherence length.
from the Eilenberger equation (36), with P x denoting the path-ordering operator. A relation similar to Eq. (53) connects the boundary matrices Q R/L with the Green's function in the far right/left region of the wire, assuming that x → ±∞. The transfer matrix satisfies certain symmetries. Along with the unitarity condition (cf. Eq. (33)) it also satisfies to the p-h symmetry: with the use of Eq. (33) one obtains and this yields We now aim to represent the Q-matrix in the scattering region in terms of the transfer matrix M and the asymptotic Eilenberger functions Q ± . We start by translating the transfer matrix relation (55) to a set of algebraic conditions 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 general produce exponentially increasing and decreasing contributions. The former are inacceptable in that they lead to exponential divergencies in the quasiclassical Green function. As is detailed in Appendix B, the requirement of a non-divergent Green function leads to the algebraic conditions while the right matrix Q R should obey the analogous relations We finally combined these equations with the transfer matrix 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 in Appendix B we obtain where M ≡ M (x R , x L ). These formulae define the starting point for an efficient numerical computation of the disordered Eilenberger function Q(x). To this end, one computes the transfer matrix M by numerical solution of corresponding system of linear first order differential equations. One next applies Eqs. (60) and (61) to obtain Q R/L . Finally, our main object of interest, Q(x), is obtained by application of M (x, x R/L ) to either Q R or Q L . So far we have considered a quantum wire connected to two superconducting terminals. However, the generalization of the method to the system shown in Fig. 1 is straightforward. The key observation is that in the limit of vanishing barrier conductance g T 1 the chiral fermion fields satisfy where φ is (energy dependent) reflection phase shift. (Here, we assumed the absence of barrier spin flip scattering, inter-channel scattering or related complications.) In the limit of asymptotically high potential barrier φ = π. Relations (62) define boundary conditions for the Eilenberger function Q L . As verified in Appendix D, these conditions assume the form of Eqs. (58), where, however, the role of Q − is taken by an "effective" matrix Q − ≡ Q − (φ) describing the tunnel junction. The explicit form of this matrix reads This matrix also satisfies Q 2 − = 1, and the relations (60,61) stay intact.
We close this section with an important statement: the Eilenberger functions Q + and Q − of the terminals define the "Majorana number" 21 of the wire as in terms of the product of two Z 2 topological invariants given by Eq. (50). It is shown in Appendix C, that for M = −1 the system supports a Majorana fermion localized in the scattering region between two terminals.

V. NUMERICS
We next turn to the discussion of our numerical results obtained for the setup shown in Fig. 1. We have solved the quasiclassical equations according to the algorithm of section IV E and from there computed the local density of states (LDoS) ν L ( ) = (2πv) −1 Re tr(σ RL z Q L ) at the left end of the wire close to the tunnel barrier. In the actual calculations we shifted the energy into the complex plane, → + iΓ. This shift accounts for the fact that in the "real" system states may escape to a continuum of lead scattering states, which gives them the status of quantum resonances of finite lifetime ∼ Γ −1 . The corresponding decay rate 22 is given by the standard golden rule expression, Γ ∼ g T δ where δ ∼ πv/2N L is the mean level spacing in the scattering region of size L. In principle, one might numerically compute the broadening by an extension of the numerical setup so as to include and extended normal metallic scattering region to the left of the tunnel barrier. This, however, would slow the performance of the numerics which is why we prefer to introduce the broadening "by hand". At any rate, the Majorana peak, present for M = −1, will acquire a finite width ∼ Γ. The DoS profiles obtained as a result of the procedure outlined above are presented and discussed in section II above.

VI. CONCLUSIONS
In this paper we have adapted the Eilenberger quasiclassical approach to the specific conditions of a class D topological quantum wire. The most striking feature of the quasiclassical approach is that the Green functions of the system are described by ordinary, and hence easily solvable differential equations. A numerical solution of these equations for given realizations of the disorder produces accurate information on the spectral properties of reasonably realistic model systems hosting Majorana boundary states. Our results confirm the predictions made in Ref. 6 for a more idealized system, viz. that disorder generates spectral peaks which can be confused for genuine Majorana peaks. It thus seems that an unambiguous detection of the Majorana might call for a measurement scheme beyond direct tunneling spectroscopy.

VII. ACKNOWLEDGMENTS
We thank C. W. J. Beenakker, F. von Oppen and M. Feigel'man for fruitful discussions. This work was supported by the collaborative research grant SFB/TR12 of the Deutsche Forschungsgemeinschaft.
We are now in position to derive Eqs. (60) and (61). To this end, we take Eqs. (58) and multiply it by the transfer matrix M from the left and by M −1 from the right. Bearing in mind Eq. (55), one obtains Adding this relation to Eqs. (59), we arrive at which in one more step yields the result (60). Eq. (61) is proven analogously.

Appendix C: Majorana mode
In this Appendix we discuss the "Majorana number" M and show how the Majorana state emerges from the quasiclassical Eilenberger function.
We consider the spectrum {E j } of Andreev bound states which follows from the poles of the Qfunction (60,61). If we denote then the energies E j are solutions of the secular equation det D(E j ) = 0. The Majorana state, if it exists, corresponds to E 0 = 0. To proceed, we introduce matricesQ ± ( ) = Q ± ( )σ RL z , and the secular matrixD( ) = D( )σ RL z . Since in the subgap interval of energies there is no distinction between the retarded and advanced Green's function, the unitarity relation (G R/A ) † = G A/R , the basic definitions Eqs. (C2) It satisfiesD T ( ) = −D(− ), and thereby guarantees that Andreev bound states appear in pairs ±E j . One may ask now whether a zero energy solution of the secular equation exists or not. At = 0, the particle hole symmetry puts tighter restrictions on the matrices Q ± , M andD (we omit the energy argument for brevity). One getsQ We thus observe that bothQ ± andD are real antisymmetric matrices of size 8N × 8N , which enables us to rewrite the secular equation Let us denote by Ω L 0 the left null space of the matrixD, i.e. any bra φ| ∈ Ω L 0 by definition satisfies φ|D = 0. SinceD is antisymmetric, the dimension of its null space is even, dim Ω L 0 = 2N . At this stage we make use of a mathematical lemma proven in Appendix A of Ref. 24: the parity of the number N can be expressed as (C7) The particle-hole symmetry (C4) implies that Det M = ±1 at zero energy. Now, the transfer matrix M (x, x ) is a continuous function of its arguments, with initial value M (x = x ) = 1. We thus conclude that Det M = +1, so that the parity of N is determined by the terminal configurations, This parity is equal to the "Majorana number" M introduced in section IV E.
Let us now focus on the most interesting case N = 1, corresponding, as we will see, to a single Majorana mode 25 . In this case we have dim Ω L 0 = 2, and the null space is spanned by two linearly independent vectors. Let φ 1 | ∈ Ω L 0 be the first basis vector. It is easy to check that φ 2 | = φ 1 |Q + ∈ Ω L 0 , and we may choose this state for the second basis vector in Ω L 0 . Indeed, for any φ 1 | ∈ Ω L 0 one has φ 1 |D = φ 1 |D = 0. Using the definition of D, Eq. (C1), we deduce that This relation enables us to evaluate φ 2 |D as and hence we proved that φ 2 | ∈ Ω L 0 . Since D is a real antisymmetric matrix, its right null space is obtained as Ω R 0 = (Ω L 0 ) † . In other words, the two kets |φ 1,2 = ( φ 1,2 )|) † satisfy the relationD|φ 1,2 = 0.
Let us look at the Eilenberger function Q R ( ) around its pole at = 0. According to Eq. (60), it can be represented by two equivalent equations, σ RL 3QR σ RL 3 = −σ RL 3 + 2(1 − Q T + )D −1 . (C14) At → 0 the inverse operator from the secular matrix has the pole structure,D −1 ∼ R/ , where the matrix R is its residue at zero energy.
At this stage it is advantageous to introduce a new bra and ket basis in the null spaces Ω L 0 and Ω R 0 , resp. The bra basis χ ± | is not orthogonal and we can denote by |η ± the ket basis, which is dual to it, χ σ |η σ = δ σσ . Using these definitions, we may formulate resolutions of unity, where σ = ± is summed over. Now let us look at the singular part of the matrixQ R , We note that the matrix P − = 1 2 (1 − Q + ) acts as the projector in the space Ω L 0 , since Similar relations hold in the ket space Ω R 0 : Equivalently, we may write Using these properties, it follows σ RL When deriving the second relation, we have taken into account that σ RL zr σ RL z = −r. The above conditions are valid for any x > x L . Sending now x → x L + 0 and using the relations (B1) written for the matrix Q L , one obtains (1 −R)(1 − Q L ) = 0, We see that these boundary conditions are equivalent to the algebraic conditions (58) if one identifies Q − = −R.
The normalizationR 2 = 1 guaranties that Q − belongs to the manifold of the quasiclassical Eilenberger functions. Transforming matrix Q − into the Majorana representation (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 characterizes the tunnel junction. In our numerical simulations we have used φ = π for all bands, which corresponds to the infinitely high barrier as compared to the energies of Andreev bound states.