Quasiclassical theory of disordered multi-channel Majorana quantum wires

Semiconductor quantum wires proximity coupled to a conventional superconductor and subjected 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–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) that resembles the Majorana peak in practically all relevant aspects. Specifically, it (i) is rigidly locked to zero energy, (ii) is of narrow width of $\mathcal {O}(\delta )$, where δ is the single-particle level spacing, (iii) carries integrated spectral weight $\mathcal {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 do show strong sensitivity to temperature, which may be due to either an intrinsic sensitivity of the peak or 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 [6], as well as those of [7, 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 this paper, we will explore the role of disorder scattering within a model much closer to the experimental setup. The price to be paid for this more realistic description is that a fully analytical treatment is out of the question. Instead, we will employ a semi-analytical approach based on the formalism of quasiclassical Green functions. Introduced in the late 1960s [9, 10], the latter has become an indispensable tool in mesoscopic superconductivity [11–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-one-dimensional topological superconductors). Specifically, we will show that in the problem at hand, the quasiclassical 'approximation' is actually very mild. Further, the quasiclassical Green 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 [6].

The rest of the paper is organized as follows. In section 2, we discuss the principal role played by disorder in the system, the idea of the quasiclassical approach and its main results. In section 3, we specify our model system, the quasi-classical approach is introduced in section 4, and in section 5 we discuss the numerical solution of the quasiclassical equations. We conclude in section 6. A number of technical details are relegated to the appendices.

A schematic view of the device currently under experimental investigation is shown in figure 1. A semiconductor quantum wire subjected to strong spin–orbit interaction is brought in contact with 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 (DoS) 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.

Zoom In Zoom Out Reset image size

In a manner to be discussed in more detail below, one contribution to the band center DoS is provided by a Majorana bound state localized at the tunnel barrier. The second, spectral contribution is generated by a conspiracy between all the 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 the 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 the order of the mean level spacing, the zero-energy 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 [8] rests on the same principal mechanism of midgap quantum interference.

Below, we will explore the phenomenon of spectral peak formation in a setting modeled to closely mimic '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 with the Fermi wavelength, yet smaller than the relevant coherence length, and this makes them suitable for quasiclassical treatment. The quasiclassical theory of the present paper is formulated in terms of the Eilenberger function $\skew3\tilde{Q} (x;\epsilon )$, 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 particle–hole degrees of freedom. It turns 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, ${\mathrm { Pf}}(\skew3\tilde{Q})$, will be seen to assume two values (±1), which locally (in space) identify the $\mathbbm{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_{\mathrm { c}}\equiv \sqrt {\Delta ^2+\mu ^2}$, where μ and Δ are chemical potential and bulk order parameter, the invariant is trivial or non-trivial, respectively. In the latter case, the wire supports a Majorana state at the barrier; in the former it does not.

Figure 2 (left) 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 curves correspond to a situation without and with the 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 figure 2 (left).

Zoom In Zoom Out Reset image size

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 the experiment,

$$\begin{equation} \frac{\mathrm{d}I}{\mathrm{d}V} = \frac{e^2\, g_T }{16\pi\hbar} \int_{-\infty}^{\infty} \frac{\partial f_{\rm F}}{\partial \epsilon}(\epsilon-V)\,\frac{\nu_{\mathrm{L}}(\epsilon)}{\nu} \mathrm{d}\epsilon, \end{equation} \tag{ 1 }$$

is essentially determined by the LDoS², i.e. the structures shown in figure 2 are expected to reflect directly in the measured signal. (Here, f_F is the Fermi distribution, ν_L is the DoS at the left barrier, ν is the DoS in the single chiral channel per unit length and V is the applied voltage.)

The profiles of the curves shown in figure 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 figure 2 (left) 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 figure 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 a dot is given by δ ≃ πv/2NL. In the presence of a magnetic field, the quasiparticle gap in the superconducting region of the wire reads $\epsilon _- = |B-\sqrt {\Delta ^2 + \mu ^2}|$. The number of subgap Andreev levels forming the effective dot is thus given by N_levels ≃ 2₋/δ.

The profiles shown in figure 2 (left) are ensemble averages, 〈ν_L〉, of the LDoS, ν_L, where the sampling was over ∼500 randomly chosen impurity configurations. To demonstrate the weakness of fluctuations, figure 2 (right) compares the average LDoS (solid line) to the 'typical' LDoS, i.e. the average $\sqrt \langle \nu _{\mathrm {L}}^2\rangle$. The relatively minor deviation between average and typical DoS demonstrates that the standard deviation

$$\begin{equation} \delta \nu(\epsilon) = \langle (\nu(\epsilon)- \langle\nu(\epsilon)\rangle)^2\rangle^{1/2} \end{equation} \tag{ 2 }$$

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 figure 3, where two unaveraged DoS profile individual disorder configurations are shown. The data are 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.

Zoom In Zoom Out Reset image size

We finally note that both the Majorana peak and the D spectral peak crucially rely on the presence and the direction 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 Majorana resonance. The class D spectral peak will also disappear if the magnetic field is rotated perpendicular to the wire (along the $\hat z$-direction). As we discussed above, under the assumption of a thin wire, L_z ≪ l_s.o., the spin–orbit term $\propto \hat s_x \partial _z$ in the Hamiltonian (4) can be neglected. Then the Bogoliubov–de-Gennes (BdG) Hamiltonian  ${\cal H}$ is invariant under rotation about the z-axis, $[{\cal H}, \hat s_z] =0$, and thus is split into two sectors, each of them belonging to symmetry class A. The wire will be in the gapped phase provided B < Δ. However, the class A one-dimensional insulator does not host boundary Majorana states, nor does the random matrix theory model of class A display any spectral anomaly at zero energy in the average DoS. Therefore, the disappearance of the conductance peak under the rotation of a magnetic field cannot be used as unambiguous evidence in favor of the Majorana particle.

In the rest of the paper, we discuss how the results summarized here were obtained by a combination of quasiclassical and numerical methods.

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 BdG Hamiltonian describing this system in the presence of disorder (section 3.1). 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 3.2), and finally transform the Hamiltonian into the Majorana basis (section 3.3). Experts in the description of topological quantum wires may proceed directly to the end of the section.

3.1. Bogoliubov–de-Gennes Hamiltonian

Introducing a four-component spinor $\Psi = (\psi _\uparrow ,\psi _\downarrow , \bar \psi _\uparrow ,\bar \psi _\downarrow )$ in the product of spin and particle–hole spaces, the BdG Hamiltonian $\cal H$ describing the system in the xz-plane reads

$$\begin{equation} {\cal H} = \frac{1}{2}\int \bar\Psi(x,z)\left( \begin{array}{@{}cc@{}} \skew3\hat{H}_0 + \hat W & i \hat s_y \Delta^* \\ -\mathrm{i} \hat s_y \Delta & -\skew3\hat{H}_0^{\mathrm{T}} - \hat W^{\mathrm{T}} \end{array}\right)\Psi(x,z)\, \mathrm{d}x\, \mathrm{d}z \end{equation} \tag{ 3 }$$

with

$$\begin{equation} \skew3\hat{H}_0 = -(\partial_x^2+\partial_z^2)/2m - \mu(x) + B(x) \hat s_x - \mathrm{i} \alpha ( \hat s_z \partial_x - \hat s_x\partial_z ). \end{equation} \tag{ 4 }$$

Here, Δ = Δ(x) is the proximity amplitude induced by the superconductor, $B\equiv \frac 12\,g \mu _{\mathrm {B}} H$, where H is the external field, we have taken into account the transverse momentum (−i∂_z) in the Rashba term, and $\hat W(x)$ is the random disorder Hamiltonian. The Pauli matrices $\hat s$ 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

$$\begin{equation} \psi_\sigma = \sum_{n\sigma} \Phi_n^\sigma(z) \psi_\sigma^{(n)} (x), \quad \bar\psi_\sigma = \sum_{n\sigma} \Phi_n^\sigma(z)^* \bar\psi_\sigma^{(n)} (x). \end{equation} \tag{ 5 }$$

Defining

$$\begin{eqnarray} &&\Psi^{(n)}=(\psi^{(n)}_{\uparrow}, \psi^{(n)}_{\downarrow}, \bar\psi^{(n)}_{\uparrow}, \bar\psi^{(n)}_{\downarrow})^{\mathrm{T}}, \quad \bar\Psi^{(n)}= (\bar\psi^{(n)}_{\uparrow}, \bar\psi^{(n)}_{\downarrow}, \psi^{(n)}_{\uparrow}, \psi^{(n)}_{\downarrow}), \end{eqnarray} \tag{ 6 }$$

the Hamiltonian then takes the form

$$\begin{eqnarray} \fl {\cal H}\!=\!\frac{1}{2}\int\!\!\bar\Psi^{(n)}\! \left( \begin{array}{@{}cc@{}} \skew3\hat{H}^{(n)}_0 \delta_{nm} + (\mathrm{i}\alpha \hat s_x\partial_z)_{nm} + \hat W_{nm} & \mathrm{i} \hat s_y \Delta^*\delta_{nm} \\ -\mathrm{i} \hat s_y \Delta\delta_{nm} & \!\!\! -(\skew3\hat{H}_0^{(n)} )^{\mathrm{T}}\delta_{nm} - (\mathrm{i}\alpha \hat s_x\partial_z)_{mn} - \hat W^{\mathrm{T}}_{mn} \end{array}\right)\Psi^{(m)}\,\mathrm{d}x, \end{eqnarray} \tag{ 7 }$$

where the one-dimensional Hamiltonian $\skew3\hat{H}^{(n)}_0$ acts in the nth band as

$$\begin{equation} \skew3\hat{H}^{(n)}_0 = -\partial_x^2/2m +\mu_z^{(n)} - \mu(x) + B(x) \hat s_x - \mathrm{i} \alpha \hat s_z \partial_x. \end{equation} \tag{ 8 }$$

In the case of an ideal waveguide the transverse wavefunctions are $\Phi _n^\sigma (z) = \sqrt {2/L_z} \sin ( n \pi z/L_z)$, so that μ⁽ⁿ⁾_z = π²n²/(2mL²_z) and the matrix elements of the spin–orbit interaction read

$$\begin{equation} h_{nm}^{\rm s.o.} = (\mathrm{i}\alpha \hat s_x\partial_z)_{nm} = -\frac{2\mathrm{i}\alpha}{L_z}\frac{nm}{n^2-m^2}\,(1-(-1)^{n+m}) \, \hat s_x. \end{equation} \tag{ 9 }$$

Let us now assume a thin wire, L_z ≲ l_so = ℏ/(mα), where l_so is the spin–orbit length³. In this case the matrix elements h^s.o._nm ≪ μ⁽ⁿ⁾_z can be treated as perturbations. To this end, we introduce a unitary transformation $\hat {\cal U}$, which brings the high-energy part of the Hamiltonian, (μ⁽ⁿ⁾_zδ_nm + h^s.o._nm), to diagonal form. Due to the weakness of the perturbation, the transformation is close to unity, $\hat {\cal U}=\exp (\mathrm { i} \hat{X} ) \simeq 1 + \mathrm {i} \hat X$, with generators $\hat X_{nm} = {\cal O}(L_z/l_{\mathrm {so}})$. To first-order perturbation theory their explicit form reads

$$\begin{equation} \hat X_{n\neq m} \simeq \frac{\mathrm{i} h^{\rm s.o.}_{nm}}{ \mu_z^{(m)} - \mu_z^{(n)}}, \quad \hat X_{nn} = 0. \end{equation} \tag{ 10 }$$

The transformation $\Psi \to \hat {\cal U} \Psi$ and $\bar \Psi \to \bar \Psi \hat {\cal U}^\dagger$ generates the approximate Hamiltonian

$$\begin{eqnarray} \fl {\cal H} \simeq \frac{1}{2}\int \bar\Psi^{(n)}(x)\left( \begin{array}{@{}cc@{}} \skew3\hat{H}^{(n)}_0 \delta_{nm} + \hat V_{nm} + \hat W_{nm} & \mathrm{i} \hat s_y \Delta\delta_{nm} + (\delta\hat\Delta)_{nm} \\ -i \hat s_y \Delta\delta_{nm} - (\delta\hat\Delta)^*_{nm} & -(\skew3\hat{H}_0^{(n)} )^{\mathrm{T}}\delta_{nm} - \hat V^{\mathrm{T}}_{mn} - \hat W^{\mathrm{T}}_{mn} \end{array}\right)\Psi^{(m)}(x)\,\mathrm{d}x. \end{eqnarray} \tag{ 11 }$$

Here we have redefined our notation for the disorder potential $\hat W\to \hat {\cal U}\hat W \hat {\cal U}^\dagger$, and introduced a small correction to the quasiparticle Hamiltonian, $\hat V = \alpha \, [\hat{X} , \hat s_z]\, \partial _x$, and to the order parameter, $\delta \hat \Delta = - \Delta [ \hat X, \hat s_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 primarily 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 $\hat W {>} {\mathrm { max}}\{\hat V,\delta \hat \Delta \}$ masks the off-diagonal matrix elements of the deterministic Hamiltonian. Having this in mind we thus omit both $\hat{V}$ and $\delta \hat \Delta$ 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 Nth 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 B² > B²_c = Δ² + (μ − μ^(N)_z)². We also assume a hierarchy of energy scales μ^(N)_z − μ^(N−1)_z ≫ E_so ≫ Δ ∼ B, where E_so = mα²/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 B² ≪ Δ² + (μ − μ⁽ⁿ⁾_z)². In the current experiments E_so ≳ Δ, i.e. our theory is on the border of applicability.

3.2. One-dimensional chiral fermions

For strong spin–orbit coupling, E_so ≫ B, each band is characterized by two Fermi momenta,

$$\begin{equation} k^{a/b}_{n} = \mp m\alpha + \sqrt{2m \left(\mu_z^{(N)}- \mu_z^{(n)} + E_{\rm so}\right)}, \end{equation} \tag{ 12 }$$

as shown in figure 4. We aim to construct a low-energy Hamiltonian describing the system at energy scales ≲E_so. To this end we represent the fields ψ⁽ⁿ⁾_σ as

$$\begin{equation} \psi_{\uparrow}^{(n)} (x) \simeq R_{a\uparrow}^{(n)}(x)\,{\rm e}^{{\rm i} k_{n}^a x} + L_{b\uparrow}^{(n)}(x)\,{\rm e}^{-{\rm i} k_{n}^b x}, \end{equation} \tag{ 13 }$$

$$\begin{equation} \psi_{\downarrow}^{(n)} (x) \simeq L_{a\downarrow}^{(n)}(x)\,{\rm e}^{-{\rm i} k_{n}^a x} + R_{b\downarrow}^{(n)}(x)\,{\rm e}^{{\rm i} k_{n}^b x} \end{equation} \tag{ 14 }$$

in terms of a superposition of right (R) and left (L) chiral fermions, and then linearize the Hamiltonian (11) 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 L-movers have spin down, while for 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.

Zoom In Zoom Out Reset image size

We next collect all R and L fields into two spinors,

$$\begin{eqnarray} &&\Psi = (R_{a\uparrow}, R_{b\downarrow}, L_{a\downarrow}, L_{b\uparrow}, \skew3\bar{R}_{a\uparrow}, \skew3\bar{R}_{b\downarrow}, \bar L_{a\downarrow}, \bar L_{b\uparrow}), \end{eqnarray} \tag{ 15 }$$

and $\bar \Psi = (\sigma _x^{\mathrm { ph}}\Psi )^{\mathrm {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

$$\begin{equation} \fl {\cal H} = \frac{1}{2}\int \bar\Psi^{(n)}(x)\left( \begin{array}{@{}cc@{}} \skew3\hat{H}_0^{(n)} \delta_{nm} + \hat W_{nm} & \mathrm{i} \Delta\,\sigma_y^{\mathrm{RL}}\otimes\sigma^{ab}_z \,\delta_{nm} \\ -\mathrm{i} \Delta\, \sigma_y^{\mathrm{RL}}\otimes\sigma^{ab}_z \,\delta_{nm} & -(\skew3\hat{H}_0^{(n)} )^{\mathrm{T}}\delta_{nm} - \hat W^{\mathrm{T}}_{mn} \end{array}\right)_{\rm ph}\Psi^{(m)}(x)\,\mathrm{d}x, \end{equation} \tag{ 16 }$$

where

$$\begin{equation} \skew3\hat{H}_0^{(n)} = \left( \begin{array}{@{}cc@{}} -\mathrm{i} v_n \partial_x - \mu & (B/2)\, ( \mathbbm{1}^{ab} + \sigma^{ab}_z)\, \delta_{Nn} \\ (B/2)\, ( \mathbbm{1}^{ab} + \sigma^{ab}_z) \, \delta_{Nn} & \mathrm{i} v_n \partial_x - \mu \end{array}\right)_{RL}. \end{equation} \tag{ 17 }$$

Here, v_n = [2m(μ^(N)_z − μ⁽ⁿ⁾_z + E_so)]^1/2/m is the Fermi velocity of the nth channel, the chemical potential μ is now defined relative to the energy μ^(N)_z, and $\hat{W} _{nm}(x)$ is a random 8 × 8 matrix satisfying the Hermiticity condition $(\hat W_{nm})^\dagger = \hat W_{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. In the case of a topological channel with n = N, the Fermi momentum k_a = 0; hence one has to keep the matrix element of magnetic field B between the R and L states. The quasiclassical approximation (i.e. the linearization of kinetic energy around the chemical potential μ_N) is, nevertheless, plausible here due to our assumption E_so ≫ max{B,Δ}. We have also used the fact that the order parameter Δ(x) can be chosen to be real.

3.3. Majorana representation

The 8N × 8N first-quantized matrix defining the Hamiltonian (16) satisfies the p–h symmetry

$$\begin{equation} \skew3\hat{H} = -\sigma_x^{\rm ph} \skew3\hat{H}^{\mathrm{T}} \sigma_x^{\rm ph}, \end{equation} \tag{ 18 }$$

which is the defining condition for a class D Hamiltonian (here, the transposition acts on the kinetic term as ∂^T_x = −∂_x). However, all that follows will be more conveniently 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

$$\begin{eqnarray*} &&\begin{array}{l} \xi_{a}^{\mathrm{R}} = (R_{a\uparrow} + \skew3\bar{R}_{a\uparrow})/\sqrt{2}, \quad \eta_{a}^{\mathrm{R}} = (R_{a\uparrow} - \skew3\bar{R}_{a\uparrow})/\sqrt{2}\mathrm{i}, \\ \xi_{b}^{\mathrm{R}} = (R_{b\downarrow} + \skew3\bar{R}_{b\downarrow})/\sqrt{2}, \quad \eta_{b}^{\mathrm{R}} = (R_{b\downarrow} - \skew3\bar{R}_{b\downarrow})/\sqrt{2}\mathrm{i} \end{array} \end{eqnarray*}$$

(and analogous relations for the L-movers, with spin reversed), which we combine into the 8N-spinor

$$\begin{equation} \tilde\chi = (\xi_{a}^{\mathrm{R}}, \xi_{b}^{\mathrm{R}}, \xi_{a}^{\mathrm{L}},\xi_{b}^{\mathrm{L}}, \eta_{a}^{\mathrm{R}}, \eta_{b}^{\mathrm{R}}, \eta_{a}^{\mathrm{L}},\eta_{b}^{\mathrm{L}})^{\mathrm{T}}. \end{equation} \tag{ 19 }$$

The two spinors χ and Ψ are related by the unitary transformation

$$\begin{eqnarray} &&\begin{array}{l} \tilde\chi = U \Psi, \quad \tilde\chi^{\mathrm{T}} = \bar\Psi U^\dagger, \\ U = \mathbbm{1}^{ab}\otimes \mathbbm{1}^{\mathrm{RL}}\otimes \displaystyle\frac{1}{\sqrt{2}}\left( \begin{array}{@{}cc@{}} 1 & 1 \\ \mathrm{i} & -\mathrm{i} \end{array}\right)_{\rm ph}. \end{array} \end{eqnarray} \tag{ 20 }$$

We combine the Majorana fields of the a and b channels as ξ^R = (ξ^R_a,ξ^R_b)^T, and reorder the spinor components as

$$\begin{equation} \chi = (\xi^{\mathrm{R}}, \eta^{\mathrm{L}}, \eta^{\mathrm{R}}, \xi^{\mathrm{L}})^{\mathrm{T}}. \end{equation} \tag{ 21 }$$

In this representation (which will be used throughout the rest of the paper), the Hamiltonian (16) takes the form

$$\begin{eqnarray} {\cal H}&=& \sum_{n,m=1}^N\frac{1}{2}\int \chi^{\mathrm{T}}_{n}(x)\,\tilde{\mathcal{H}}_{nm}(x)\,\chi_{m}(x)\,\mathrm{d}x,\nonumber\\ \tilde{\mathcal{H}}_{nm}(x)&\equiv& \left( \begin{array}{@{}cc@{}} \skew3\hat{H}_-^{(n)} \delta_{nm} +\mathrm{i} \hat W_{nm}^{--} & \mathrm{i} \mu\,\sigma_z^{\mathrm{RL}}\otimes\mathbbm{1}^{ab}\, \delta_{nm} +\mathrm{ i} \hat W_{nm}^{-+} \\ -\mathrm{ i} \mu\,\sigma_z^{\mathrm{RL}}\otimes\mathbbm{1}^{ab}\, \delta_{nm} + \mathrm{i} \hat W_{nm}^{+-} & \skew3\hat{H}_+^{(n)} \delta_{nm} + \mathrm{i} \hat W_{nm}^{++} \end{array}\right), \end{eqnarray} \tag{ 22 }$$

where the deterministic part reads

$$\begin{equation} \skew3\hat{H}^{(n)}_\pm = -\mathrm{i} v_n\, \sigma_z^{\mathrm{RL}} \otimes\mathbbm{1}^{ab}\,\partial_x - \sigma_y^{\mathrm{RL}}\otimes\hat\Delta_\pm^{(n)}, \end{equation} \tag{ 23 }$$

$$\begin{equation} \Delta_\pm^{(n)} = \Delta\,\sigma_z^{ab} \pm (B/2)( \mathbbm{1}^{ab} + \sigma^{ab}_z) \, \delta_{Nn} \end{equation} \tag{ 24 }$$

and the random matrices are constructed as

$$\begin{equation} \hat W^{--} = \left( \begin{array}{@{}cc@{}} \hat w_2^{\mathrm{RR}} & - \hat w_1^{\mathrm{RL}} \\ \hat w_1^{\mathrm{LR}} & \hat w_2^{\mathrm{LL}} \end{array} \right)_{ab}, \quad \hat W^{-+} = \left( \begin{array}{@{}cc@{}} -\hat w_1^{\mathrm{RR}} & \hat w_2^{\mathrm{RL}} \\ \hat w_2^{\mathrm{LR}} & \hat w_1^{\mathrm{LL}} \end{array} \right)_{ab} \end{equation} \tag{ 25 }$$

in terms of real (symmetric) and imaginary (antisymmetric) parts of the random matrices $(\hat W)^{\mathrm {RL}} = \hat{w} _1^{\mathrm {RL}} + \mathrm {i}\, \hat w_2^{\mathrm {RL}}$, etc. The remaining blocks of the $\hat W$-matrix are defined by

$$\begin{equation} \hat W^{+-} = -\sigma_z^{ab}\,\hat W^{-+}\,\sigma_z^{ab}, \quad \hat W^{++} = \sigma_z^{ab}\,\hat W^{--}\,\sigma_z^{ab}.\end{equation}\noindent \tag{ 26 }$$

In the Majorana basis (21), the class D symmetry is expressed through the antisymmetry ${\skew3\hat{H}}^{\mathrm {T}}=-{\skew3\hat{H}}$.

We finally specify the statistics of disorder. We choose $\hat W(x)$ to be a δ-correlated and Gaussian distributed random matrix of size 8N × 8N with a zero mean value $\langle \hat W(x) \rangle$ and variance

$$\begin{eqnarray} &&\begin{array}{l} \langle w^{ij}_1(x) w^{i'j'}_1(x')\rangle =\displaystyle\frac{\gamma_{{\rm w}}}{2}\delta(x-x') (\delta_{ii'}\delta_{j'j'} + \delta_{ij'}\delta_{ji'}), \\ \langle w^{ij}_2(x) w^{i'j'}_2(x')\rangle =\displaystyle\frac{\gamma_{{\rm w}}}{2}\delta(x-x') (\delta_{ii'}\delta_{jj'} - \delta_{ij'}\delta_{ji'}). \end{array} \end{eqnarray} \tag{ 27 }$$

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 time reversal and spin rotation symmetry (e.g. random spin-flip scattering caused by random spin–orbit terms is included in (27)). The strength of disorder set by the coefficient γ_w translates into the 'golden rule' scattering rate

$$\begin{equation} \tau^{-1} = 2\gamma_{{\rm w}}^2 \sum_{n=1}^N (1/v_n) \end{equation} \tag{ 28 }$$

of the normal conducting (i.e. superconductor decoupled) quantum wire.

The kinetic term in the low-energy Hamiltonian (22) that 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  [19, 20]. After introducing the basics of the method (section 4.1), we construct the Eilenberger Q-function in the limit of a single clean 'topological channel' (section 4.2) and discuss the resulting DoS (section 4.3). In section 4.4, we define the $\mathbbm{Z}_2$ topological invariant in terms of the Q-matrix. Section 4.5 outlines the general construction of the solution to the Eilenberger equation in the inhomogeneous disordered wire with boundary conditions.

4.1. Eilenberger method

We start by defining

$$\begin{equation} g^{\rm R/A}_{\epsilon,nm}(x,x') = \sqrt{v_n}\, G^{\rm R/A}_{\epsilon,nm}(x,x')\, \sigma_z^{\mathrm{RL}} \, \sqrt{v_m}, \end{equation} \tag{ 29 }$$

where $G_\epsilon ^{\mathrm { R/A}}(x,x') \equiv \langle x | (\epsilon \pm \mathrm {i}0 - \skew3\tilde{H})^{-1}| x' \rangle$ are the retarded and advanced Green functions of the system. Under transposition (which in our current representation represents the particle–hole symmetry), the function g behaves as

$$\begin{equation} g_\epsilon^{\mathrm{R}}(x,x') = -\sigma_z^{\mathrm{RL}}[ g_{-\epsilon}^{\mathrm{A}}(x',x)]^{\mathrm{T}} \sigma_z^{\mathrm{RL}}, \end{equation} \tag{ 30 }$$

i.e. advanced and retarded Green functions get interchanged. It is not hard to derive the two mutually adjoint differential equations

$$\begin{eqnarray} &&\begin{array}{l} \partial_{x}\, g_\epsilon(x, x') +\mathcal{L}_\epsilon g_\epsilon(x, x') = -\mathrm{i} \delta(x-x'),\\ \partial_{x'} \, g_\epsilon(x, x') - g_\epsilon(x, x') \mathcal{L}_\epsilon = \mathrm{i}\delta(x-x') \end{array} \end{eqnarray} \tag{ 31 }$$

describing the dynamical evolution of g. Here,

$$\begin{equation} \mathcal{L}_\epsilon\equiv -\mathrm{i} (\hat{\omega} - {\cal P} ), \end{equation} \tag{ 32 }$$

where matrix $\hat \omega$ has elements

$$\begin{equation} (\hat{\omega})_{nm} \equiv (\epsilon/v_n)\, \sigma_z^{\mathrm{RL}} \,\,\delta_{nm} \end{equation} \tag{ 33 }$$

and the operator $\cal P$ is related to the Hamiltonian matrix $\tilde {\mathcal {H}}$ as

$$\begin{equation} \mathcal{P}_{nm}= \mathrm{i}\partial_x + v_n^{-1/2}\,(\sigma_z^{\mathrm{RL}} \skew3\tilde{H}_{nm} )\, v_m^{-1/2}. \end{equation} \tag{ 34 }$$

Due to the antisymmetry of $\tilde {\mathcal {H}}=-\tilde {\mathcal {H}}^{\mathrm {T}}$, the operator $\mathcal {P}$ obeys the particle–hole symmetry

$$\begin{equation} \mathcal{P}=-\sigma_z^{\mathrm{RL}} \mathcal{P}^{\mathrm{T}} \sigma_z^{\mathrm{RL}}. \end{equation} \tag{ 35 }$$

We next define the Eilenberger function as

$$\begin{equation} Q_\epsilon (x) = \lim_{x'\to x}[ 2\mathrm{i} \,g_\epsilon(x,x') - {\rm sgn}(x-x')], \end{equation} \tag{ 36 }$$

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 equations (31). Subtracting the two equations in (31), we then obtain the Eilenberger equation of motion

$$\begin{equation} \partial_x Q_\epsilon(x) + [\mathcal{L}_\epsilon, Q_\epsilon(x)] = 0. \end{equation} \tag{ 37 }$$

The Eilenberger function Q obeys the particle–hole symmetry

$$\begin{equation} \sigma_z^{\mathrm{RL}} \, Q_{-\epsilon}^{\mathrm{T}}(x)\, \sigma_z^{\mathrm{RL}} = - Q_\epsilon(x) \end{equation} \tag{ 38 }$$

and the normalization condition $Q_\epsilon ^2(x) = \mathbbm{1}$, where $\mathbbm{1}$ is the unit matrix (the latter condition can be checked by verifying that equation (37) preserves the normalization Q² = const.). The unit value of the normalization constant is fixed by the jump height of the sgn-function in (36). We finally note that the operator $\cal P$ can be straightforwardly constructed from equation (34); however, for our present purposes, we need not state its explicit form in generality.

4.2. Eilenberger function in the clean limit

As a warmup, we apply the quasiclassical approach to the limit ($\hat{W} = 0$) of an infinite clean quantum wire subjected 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

$$\begin{equation} L_\epsilon = -\mathrm{i} \left( \begin{array}{@{}cc@{}} \epsilon\, \sigma_z^{\mathrm{RL}} - \mathrm{i}\Delta_{-}\sigma_x^{\mathrm{RL}} & -\mathrm{i}\mu\, \mathbbm{1}^{\mathrm{RL}} \\ \mathrm{i}\mu\,\mathbbm{1}^{\mathrm{RL}} & \epsilon\, \sigma_z^{\mathrm{RL}} - \mathrm{i}\Delta_{+}\sigma_x^{\mathrm{RL}} \end{array} \right) \end{equation} \tag{ 39 }$$

as the reduction of the general operator $\mathcal {L}_\epsilon$ to a single channel. The solution Q is then determined by the relations [Q,L] = 0 and $Q_{\epsilon }^2=\mathbbm{1}$. To solve these equations, we assume L to be diagonalized as

$$\begin{equation} L_\epsilon = T (\hat \lambda\otimes \sigma_z^{\mathrm{RL}}) T^{-1}, \quad \hat \lambda \equiv \left( \begin{array}{@{}cc@{}} \lambda_+ & \\ & \lambda_- \end{array} \right)\!, \end{equation} \tag{ 40 }$$

where the exact form of the (non-unitary) transformation matrix T will not be needed and

$$\begin{eqnarray} &&\fl \lambda_\pm =\sqrt{\lambda_0\pm \lambda},\quad \lambda_0 = B^2 + \Delta^2 - \mu^2 - \epsilon^2, \quad \lambda = 2\sqrt{ B^2\Delta^2 - (\Delta^2-\epsilon^2)\mu^2} \end{eqnarray} \tag{ 41 }$$

are the eigenvalues. The defining equations for Q are then solved by matrices of the form

$$\begin{equation} Q_\epsilon = T\Lambda T^{-1}, \end{equation} \tag{ 42 }$$

where $\Lambda =(\pm 1,\ldots ,\pm 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 function (opposite for advanced) is given by

$$\begin{equation} \quad \Lambda =\sigma_z^{\mathrm{RL}}. \end{equation} \tag{ 43 }$$

A more explicit derivation of this structure is detailed in appendix A.

4.3. Clean density of states

The DoS in the bulk of the topological wire is given by ν() = (2πv)⁻¹Re tr σ^RL_zQ. 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

$$\begin{equation} Q_\epsilon = \frac{L_\epsilon}{\lambda_+}\, P^+ + \frac{L_\epsilon}{\lambda_-} \, P^-, \end{equation} \tag{ 44 }$$

where P⁺ and P⁻ are projectors on the space of L-eigenstates with eigenvalues ±λ₊ and ±λ₋, respectively:

$$\begin{equation} P^+=T\mathrm{diag}(\mathbbm{1}_2,0)T^{-1},\quad P^-=T\mathrm{diag}(0,\mathbbm{1}_2)T^{-1}. \end{equation} \tag{ 45 }$$

That this representation faithfully represents the matrix Q is checked by application of (44) in the eigenbasis where all matrices assume a diagonal form. It remains to obtain a representation of P^± that does not make explicit reference to the diagonalizing matrices T. To this end, note that (eigenrepresentation understood) $L_\epsilon ^2=\mathrm {diag}(\lambda _+^2 \mathbbm{1}_2,\lambda _-^2\mathbbm{1}_2)=\lambda _0\mathbbm{1}_4+\lambda P^+-\lambda P^-$. This equation can be straightforwardly solved as

$$\begin{equation} P^\pm={1\over 2}\left(\mathbbm{1}_4\pm {1\over \lambda} (L_\epsilon^2 - \lambda_0\mathbbm{1}_4)\right). \end{equation}\noindent \tag{ 46 }$$

Substituting this expression into equation (44), we obtain a representation of Q that makes reference only to the operator (39), and its eigenvalues. Computing the trace, we obtain DoS profiles as shown in figure 5. 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\leftrightarrow -\mathrm {i}\partial _x$) and the 'potential' ($\mathcal {L}$). On the other hand, we computed the DoS from Q as determined by $\mathcal {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=\mathbbm{1}$. Indeed, the conservation of the constraint, and the unit value of the normalization are consequences of the linearity of the derivative operator in (31), which in this way co-determines the structure of Q.

Zoom In Zoom Out Reset image size

Inspection of equation (44) shows that the DoS contains singularities at the zeros λ^±() = 0, which are located at

$$\begin{equation} \epsilon_{\pm} = \vert B\pm \sqrt{\Delta^2 + \mu^2}\vert. \end{equation} \tag{ 47 }$$

Let us assume that B > Δ. Figure 5(a) shows the ensuing DoS profile, along with the underlying dispersion relation for a value of the chemical potential μ < μ_c, where

$$\begin{equation} \mu_{\rm c} = \sqrt{B^2 - \Delta^2} \end{equation} \tag{ 48 }$$

defines a critical value where the lower of the DoS singularities, ₋, touches zero and the band gap closes (figure 5(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

$$\begin{equation} \mu^* = \sqrt{B^2 + B \sqrt{B^2 + 4\Delta^2}}/\sqrt{2}, \end{equation} \tag{ 49 }$$

the lower band ₋(k) develops an extremum at finite values of k that manifests in a third van-Hove singularity at the energy

$$\begin{equation} \epsilon_0 = \Delta\sqrt{1-B^2/\mu^2}, \end{equation} \tag{ 50 }$$

as shown in panel (d) of figure 5. This energy minimum corresponds to the Fermi momentum k_a of the 'topological' channel that effectively becomes non-zero at high chemical potential μ (see equation (12)).

4.4. $\mathbbm{Z}_2$ topological invariant

The symmetry equation (38) implies that at zero energy the product σ^RL_zQ₌₀ 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 $\sigma ^{\mathrm {RL}}_z=\sigma _z^{\mathrm {RL}}\otimes \mathbbm{1}_2$. Consequently, the Pfaffian of σ^RL_zQ₌₀—which squares to the determinant of that matrix—may take one of two values, ±1. This motivates the definition of the topological index

$$\begin{equation} N_{\rm top} = {\rm Pf} (\sigma_z^{\mathrm{RL}}\, Q_{\epsilon=0}) = \left\{ \begin{array}{@{}lc} +1, & \mu > \mu_{\rm c}, \\ -1, & \mu < \mu_{\rm c} \end{array} \right. \end{equation} \tag{ 51 }$$

distinguishing between the two phases. Computing N_top from (44), we find the index structure stated in (51). (Right at the critical point μ = μ_c, the matrix Q₌₀ becomes singular and the index cannot be defined.)

4.5. 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 figure 6. The terminals indicated at the left and right represent superconducting regions, assumed to be non-disordered for simplicity. (This is an inconsequential assumption provided the rate of disorder scattering τ⁻¹ ∼ Nγ²_w〈1/v_n〉 does not exceed the energy gaps (₋ or ₀) 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 $\mathcal {L}$, we will show how the left and right asymptotic configurations 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.

Zoom In Zoom Out Reset image size

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 (figure 6). These constants set asymptotic values of the Q-matrix,

$$\begin{equation} Q_\epsilon(x\to -\infty) \equiv Q_{-}, \quad Q_\epsilon(x\to +\infty) \equiv Q_+, \end{equation} \tag{ 52 }$$

where Q_± are constructed using the results of section 4.2 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, respectively, as

$$\begin{equation} Q_{\mathrm{R}} = Q(x_{\mathrm{R}}), \quad Q_{\mathrm{L}}=Q(x_{\mathrm{L}}), \end{equation}\noindent \tag{ 53 }$$

where x_R = −x_L = L/2. These two configurations are related by a transfer matrix,

$$\begin{equation} Q_{\mathrm{R}} = M(x_{\mathrm{R}}, x_{\mathrm{L}})\, Q_{\mathrm{L}}\, M^{-1}(x_{\mathrm{R}}, x_{\mathrm{L}}). \end{equation}\noindent \tag{ 54 }$$

For arbitrary positions x and x', the formal expression for the transfer matrix M(x,x') at a given energy follows from the Eilenberger equation (37),

$$\begin{equation} M_\epsilon(x,x') = {\cal P}_x \exp\left\{\mathrm{i}\int_{x}^{x'} \left[\hat\omega- {\cal P}(y)\right] \mathrm{d}y \right\} \end{equation} \tag{ 55 }$$

with ${\cal P}_x$ denoting the path-ordering operator. A relation similar to equation (54) connects the boundary matrices Q_R/L with the Green function in the far right/left region of the wire,

$$\begin{equation} Q_{\pm} = M(x, x_{\mathrm{R/L}})\, Q_{\mathrm{R/L}}\, M^{-1}(x,x_{\mathrm{R/L}}), \end{equation} \tag{ 56 }$$

assuming that x → ± ∞. The transfer matrix satisfies certain symmetries. Along with the unitarity condition (cf equation (34))

$$\begin{equation} \sigma_z^{\mathrm{RL}} \,M_\epsilon^\dagger\,\sigma_z^{\mathrm{RL}} = M_\epsilon^{-1}, \end{equation} \tag{ 57 }$$

it also satisfies the p–h symmetry: with the use of equation (34) one obtains

$$\begin{equation} \mathcal{L}_\epsilon^{\mathrm{T}} = -\sigma_z^{\mathrm{RL}} \mathcal{L}_{-\epsilon}\sigma_z^{\mathrm{RL}}, \end{equation} \tag{ 58 }$$

and this yields

$$\begin{equation} \sigma_z^{\mathrm{RL}}\, M_{-\epsilon}^{\mathrm{T}} \, \sigma_z^{\mathrm{RL}} = M_{\epsilon}^{-1}. \end{equation} \tag{ 59 }$$

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 (56) to a set of algebraic conditions relating the matrices Q_R/L to Q_±. To this end, note that the action of the non-unitary (cf equation (57)) 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

$$\begin{eqnarray} &&\begin{array}{l} (\mathbbm{1}+Q_-)(\mathbbm{1}-Q_{\mathrm{L}})=0, \\ (\mathbbm{1}+Q_{\mathrm{L}})(\mathbbm{1}-Q_-)=0, \end{array} \end{eqnarray} \tag{ 60 }$$

while the right matrix Q_R should obey the analogous relations

$$\begin{eqnarray} &&\begin{array}{l} (\mathbbm{1}-Q_+)(\mathbbm{1}+Q_{\mathrm{R}})=0, \\ (\mathbbm{1}-Q_{\mathrm{R}})(\mathbbm{1}+Q_+)=0. \end{array} \end{eqnarray} \tag{ 61 }$$

We finally combined these equations with the transfer matrix relation (56) 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

$$\begin{equation} Q_{\mathrm{R}} = \mathbbm{1} + \frac{2}{Q_+ + M Q_- M^{-1}}\,(\mathbbm{1}- Q_+), \end{equation} \tag{ 62 }$$

$$\begin{equation} Q_{\mathrm{L}} = \mathbbm{1} + (\mathbbm{1}- Q_-)\,\frac{2}{Q_- + M^{-1} Q_+ M}, \end{equation} \tag{ 63 }$$

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 the corresponding system of linear first-order differential equations. One next applies equations (62) and (63) 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 figure 1 is straightforward. The key observation is that in the limit of vanishing barrier conductance g_T  ≪ 1, the chiral fermion fields satisfy

$$\begin{equation} R_{a\uparrow}(x_{\mathrm{L}}) = \mathrm{e}^{\mathrm{i}\phi} L_{b\uparrow}(x_{\mathrm{L}}),\quad R_{b\downarrow}(x_{\mathrm{L}}) = \mathrm{e}^{\mathrm{i}\phi} L_{a\downarrow}(x_{\mathrm{L}}), \end{equation} \tag{ 64 }$$

where ϕ is the (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 (64) define boundary conditions for the Eilenberger function Q_L. As verified in appendix D, these conditions assume the form of equations (60), 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

$$\begin{equation} Q_- = (\sin(\phi)\sigma_z^{\rm ph} + \cos(\phi)\sigma_x^{\rm ph} )\otimes \sigma_x^{\mathrm{RL}}\otimes\sigma_x^{ab}. \end{equation} \tag{ 65 }$$

This matrix also satisfies $Q_-^2=\mathbbm{1}$, and the relations (62) and (63) 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

$$\begin{equation} {\cal M} = {\rm Pf} (\sigma_z^{\mathrm{RL}}\, Q_+) \, {\rm Pf} (\sigma_z^{\mathrm{RL}}\, Q_-) \end{equation} \tag{ 66 }$$

in terms of the product of two $\mathbbm{Z}_2$ topological invariants given by equation (51). It is shown in appendix C that for ${\cal M}=-1$ the system supports a Majorana fermion localized in the scattering region between two terminals.

We next turn to the discussion of our numerical results obtained for the setup shown in figure 1. We have solved the quasiclassical equations according to the algorithm of section 4.5 and from there computed the LDoS ν_L() = (2πv)⁻¹Re tr(σ^RL_zQ_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 ∼Γ⁻¹. The corresponding decay rate [22] is given by the standard golden rule expression Γ ∼ g_T δ, where δ ∼ πv/2NL 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 extend the normal metallic scattering region to the left of the tunnel barrier. This, however, would slow down the performance of the numerics, which is why we prefer to introduce the broadening 'by hand'. At any rate, the Majorana peak, present for ${\cal 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 2 above.

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 [6] for a more idealized system, viz. that disorder generates spectral peaks that 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.

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 no. SFB/TR12 of the Deutsche Forschungsgemeinschaft.

The goal of this appendix is to derive the structure (43) by explicit calculation. Our starting point is the momentum representation of the equations of motion for the quasiclassical Green function g(x,x';) ≡ g(x − x';),

$$\begin{equation} \left( - p + \frac{\mathrm{i}}{v} L_\epsilon \right)\,g(p,\epsilon)=1. \end{equation} \tag{ A.1 }$$

Assuming a diagonal representation as in (40), we obtain

$$\begin{equation} g(p,\epsilon) = - \,T {1\over p-\mathrm{i}{\hat \lambda\over v}\otimes \sigma_z^{\mathrm{RL}}} T^{-1}. \end{equation} \tag{ A.2 }$$

The inspection of eigenvalues λ_± shows that for the retarded Green function, i.e. at  →  + i0, one has Re (λ_±) > 0. This gives

$$\begin{eqnarray} g(x,x';\epsilon)&=T\int {\mathrm{d}p\over 2\pi}\frac{\mathrm{e}^{\mathrm{i}p x}}{p - \mathrm{i} {\hat \lambda\over v} \otimes \sigma_z^{\mathrm{RL}}} T^{-1} \\ &=-{\mathrm{i}\over 2} T (\sigma_z^{\mathrm{RL}} +\mathrm{sgn}\, (x-x')) \,\mathrm{e}^{-{\hat \lambda \over v}|x-x'|}\, T^{-1}.\nonumber \end{eqnarray} \tag{ A.3 }$$

Sending x' → x and comparing to (36), we obtain the identifications (42) and (43).

In this appendix, we derive the boundary conditions (60), (61) for the Q-matrix. The method we use is adapted from the circuit theory of [20, 23]. For simplicity, we consider the 4 × 4 matrices describing individual channels; generalization to many coupled channels is straightforward.

Consider the quasiclassical function g(x,x') ≡ g(x,x';) to the right of the right terminal, x,x' ⩾ x_R. By assumption, the generator L | _x>xR is constant in space, and this means that equations (31) defining g admit the solutions

$$\begin{eqnarray} &&\begin{array}{l} g(x,x_{\mathrm{R}}) = \mathrm{e}^{-L_\epsilon (x-x_{\mathrm{R}})} g(x_{\mathrm{R}}+0,x_{\mathrm{R}}),\\ g(x_{\mathrm{R}},x) = g(x_{\mathrm{R}},x_{\mathrm{R}}+0)\,\mathrm{e}^{L_\epsilon (x-x_{\mathrm{R}})}, \end{array} \end{eqnarray} \tag{ B.1 }$$

where we have again set v = 1 for notational simplicity. According to equation (A.3), we have

$$\begin{eqnarray} &&\begin{array}{l} g(x_{\mathrm{R}}+0,x_{\mathrm{R}}) = -\displaystyle{\mathrm{i}\over 2}(Q_{\mathrm{R}} + \mathbbm{1}), \\ g(x_{\mathrm{R}},x_{\mathrm{R}}+0) = -\displaystyle{\mathrm{i}\over 2}(Q_{\mathrm{R}} - \mathbbm{1}) \end{array} \end{eqnarray} \tag{ B.2 }$$

so that

$$\begin{eqnarray} &&\begin{array}{l} g(x,x_{\mathrm{R}}) = -\displaystyle{\mathrm{i}\over 2}\mathrm{e}^{-L_\epsilon (x-x_{\mathrm{R}})} (Q_{\mathrm{R}} + \mathbbm{1}),\\ g(x_{\mathrm{R}},x) = -\displaystyle{\mathrm{i}\over 2} (Q_{\mathrm{R}} + \mathbbm{1}) \mathrm{e}^{L_\epsilon (x-x_{\mathrm{R}})}. \end{array} \end{eqnarray} \tag{ B.3 }$$

Defining $\bar g = T^{-1} g T$ and $\skew3\bar{Q}_{\mathrm {R}} = T^{-1} Q_{\mathrm {R}} T$, we next transform to the representation (40) to arrive at

$$\begin{eqnarray} &&\begin{array}{l} \bar g(x,x_{\mathrm{R}}) = -\displaystyle{\mathrm{i}\over 2}\mathrm{e}^{-\hat \lambda \otimes \sigma_z^{\mathrm{RL}} (x-x_{\mathrm{R}})} (\skew3\bar{Q}_{\mathrm{R}} + \mathbbm{1}), \\ \bar g(x_{\mathrm{R}},x) = -\displaystyle{\mathrm{i}\over 2} (\skew3\bar{Q}_{\mathrm{R}} + \mathbbm{1}) \, \mathrm{e}^{\hat \lambda \otimes \sigma_z^{\mathrm{RL}} (x-x_{\mathrm{R}})}. \end{array} \end{eqnarray} \tag{ B.4 }$$

The key observation now is that the matrix functions $\bar g$ must remain finite as x → ∞. Due to the positivity of the matrix $\hat \lambda$, this is equivalent to the condition $(\mathbbm{1}-\sigma _z^{\mathrm {RL}})(\hat{Q} _{\mathrm {R}}+\mathbbm{1}) = (\mathbbm{1}-\sigma _z^{\mathrm {RL}}\otimes \mathbbm{1})(\hat{Q} _{\mathrm {R}}+\mathbbm{1}) =(\mathbbm{1}-\Lambda )(\hat{Q} _{\mathrm {R}}+\mathbbm{1}) =0$. By the same token, we have $(\skew3\bar{Q}_{\mathrm {R}} +\mathbbm{1})(\mathbbm{1}+\Lambda )=0$. Finally, noting that TΛT⁻¹ = Q₊ represents the asymptotic form of the Q-matrix deep in the superconductor, we arrive at equations (61). Equations (60) are shown in an analogous way.

We are now in a position to derive equations (62) and (63). To this end, we take equations (60) and multiply it by the transfer matrix M from the left and by M⁻¹ from the right. Bearing in mind equation (56), one obtains

$$\begin{equation} \mathbbm{1} - M Q_- M^{-1}\, Q_{\mathrm{R}} + M Q_- M^{-1} - Q_+ = 0. \end{equation} \tag{ B.5 }$$

Adding this relation to equations (61), we arrive at

$$\begin{equation} 2\cdot\mathbbm{1}-M Q_- M^{-1}\, Q_{\mathrm{R}} - Q_+ Q_{\mathrm{R}} + M Q_- M^{-1} - Q_+ = 0, \end{equation} \tag{ B.6 }$$

which in one more step yields the result (62). Equation (63) is proven analogously.

In this appendix, we discuss the 'Majorana number' ${\cal 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 Q-function (62), (63). If we denote

$$\begin{equation} {\cal D}(\epsilon) = Q_+(\epsilon) + M(\epsilon) Q_-(\epsilon) M^{-1}(\epsilon), \end{equation} \tag{ C.1 }$$

then the energies E_j are solutions of the secular equation ${\mathrm { det}}\,{\cal D}(E_j)=0$. The Majorana state, if it exists, corresponds to E₀ = 0.

To proceed, we introduce matrices $\skew3\tilde{Q}_{\pm }(\epsilon ) = Q_{\pm }(\epsilon )\sigma _z^{\mathrm { RL}}$, and the secular matrix $\tilde {\cal D}(\epsilon ) = {\cal D}(\epsilon ) \sigma _z^{{\mathrm { RL}}}$. Since in the subgap interval of energies there is no distinction between the retarded and advanced Green function, the unitarity relation (G^R/A)^† = G^A/R, the basic definitions equations (29) and (36) imply $\skew3\tilde{Q}_{\pm }^\dagger (\epsilon ) = -\skew3\tilde{Q}_{\pm }(\epsilon )$. Further, particle–hole symmetry (38) yields $\skew3\tilde{Q}_{\pm }^{\mathrm {T}}(-\epsilon ) = -\skew3\tilde{Q}_\pm (\epsilon )$. Taking into account the class D symmetry of the transfer matrix, equation (59), the secular matrix takes a form

$$\begin{equation} \tilde {\cal D}(\epsilon) = \skew3\tilde{Q}_+(\epsilon) + M(\epsilon) \skew3\tilde{Q}_-(\epsilon) M^{\mathrm{T}}(-\epsilon). \end{equation} \tag{ C.2 }$$

It satisfies $\tilde {\cal D}^{\mathrm {T}}(\epsilon ) = - \tilde {\cal D}(-\epsilon )$, 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 $\skew3\tilde{Q}_{\pm }$, M and $\tilde {\cal D}$ (we omit the energy argument for brevity). One obtains

$$\begin{equation} \skew3\tilde{Q}_{\pm}^{\mathrm{T}} = - \skew3\tilde{Q}_{\pm}, \quad \skew3\tilde{Q}_{\pm}^* = \skew3\tilde{Q}_{\pm}, \end{equation} \tag{ C.3 }$$

$$\begin{equation} \sigma_z^{\rm RL} M^{\mathrm{T}}\sigma_z^{\rm RL} = M^{-1}, \quad M^* = M, \end{equation} \tag{ C.4 }$$

$$\begin{equation} \tilde {\cal D}^{\mathrm{T}} = - \tilde{\cal D}, \quad \tilde {\cal D}^* = \tilde{\cal D}. \end{equation} \tag{ C.5 }$$

We thus observe that both $\skew3\tilde{Q}_{\pm }$ and $\tilde {\cal D}$ are real antisymmetric matrices of size 8N × 8N, which enables us to rewrite the secular equation in terms of a Pfaffian

$$\begin{equation} {\rm Det}\,\tilde {\cal D} = [ {\rm Pf}(\tilde {\cal D})]^2 = [ {\rm Pf}(\skew3\tilde{Q}_+ + M \skew3\tilde{Q}_- M^{\mathrm{T}})]^2 = 0. \end{equation} \tag{ C.6 }$$

Let us denote by Ω^L₀ the left null space of the matrix $\tilde {\cal D}$, i.e. any bra 〈ϕ|∈Ω^L₀ by definition satisfies $\langle \phi | \tilde {\cal D} = 0$. Since $\tilde {\cal D}$ is antisymmetric, the dimension of its null space is even, ${\mathrm { dim}}\, \Omega ^{\mathrm {L}}_0 = 2 {\cal N}$. At this stage, we make use of a mathematical lemma proven in appendix A of [24]: the parity of the number ${\cal N}$ can be expressed as

$$\begin{equation} (-1)^{\cal N} = {\rm Pf}\,\skew3\tilde{Q}_+\, {\rm Pf}\,[M \, \skew3\tilde{Q}_-\, M^{\mathrm{T}}] = {\rm Det}\, M\, {\rm Pf}\,\skew3\tilde{Q}_-\, {\rm Pf}\,\skew3\tilde{Q}_+. \end{equation} \tag{ C.7 }$$

The particle–hole symmetry (C.4) 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') =\mathbbm{1}$. We thus conclude that Det M = + 1, so that the parity of ${\cal N}$ is determined by the terminal configurations,

$$\begin{equation} (-1)^{\cal N} = {\rm Pf}\,\skew3\tilde{Q}_-\, {\rm Pf}\, \skew3\tilde{Q}_+. \end{equation} \tag{ C.8 }$$

This parity is equal to the 'Majorana number' ${\cal M}$ introduced in section 4.5.

Let us now focus on the most interesting case ${\cal N}=1$, corresponding, as we will see, to a single Majorana mode⁴. In this case, we have dim Ω^L₀ = 2, and the null space is spanned by two linearly independent vectors. Let 〈ϕ₁|∈Ω^L₀ be the first basis vector. It is easy to check that 〈ϕ₂| = 〈ϕ₁|Q₊∈Ω^L₀, and we may choose this state for the second basis vector in Ω^L₀. Indeed, for any 〈ϕ₁|∈Ω^L₀ one has $\langle \phi _1| \tilde {\cal D} = \langle \phi _1| {\cal D} = 0$. Using the definition of ${\cal D}$, equation (C.1), we deduce that

$$\begin{equation} \langle \phi_1| Q_+ M = - \langle \phi_1| M Q_-. \end{equation} \tag{ C.9 }$$

This relation enables us to evaluate $\langle \phi _2| {\cal D}$ as

$$\begin{eqnarray} \langle \phi_2| {\cal D} &=& \langle \phi_1| \left( \mathbbm{1} + Q_+ M Q_- M^{-1}\right) = \langle \phi_1| + \left( \langle \phi_1| Q_+ M \right)\, Q_- M^{-1} \nonumber \\ &=& \langle \phi_1|- \langle \phi_1| M Q_- Q_- M^{-1} = 0 \end{eqnarray} \tag{ C.10 }$$

and hence we proved that 〈ϕ₂|∈Ω^L₀. Since ${\cal D}$ is a real antisymmetric matrix, its right null space is obtained as Ω^R₀ = (Ω^L₀)^†. In other words, the two kets |ϕ_1,2〉 = (〈ϕ_1,2)|)^† satisfy the relation $\tilde {\cal D} |\phi _{1,2}\rangle = 0$.

Let us look at the Eilenberger function Q_R() around its pole at  = 0. According to equation (62), it can be represented by two equivalent equations,

$$\begin{equation} Q_{\mathrm{R}} = \mathbbm{1} + 2{\cal D}^{-1} (\mathbbm{1} -Q_+), \end{equation} \tag{ C.11 }$$

$$\begin{equation} Q_{\mathrm{R}} = -\mathbbm{1} + 2(\mathbbm{1} + Q_+){\cal D}^{-1}. \end{equation} \tag{ C.12 }$$

Multiplication by σ^RL_z yields

$$\begin{equation} \sigma_z^{\rm RL} \skew3\tilde{Q}_{\mathrm{R}} \sigma_z^{\rm RL} = \sigma_z^{\rm RL}+ 2\tilde {\cal D}^{-1} (\mathbbm{1} -Q_+), \end{equation} \tag{ C.13 }$$

$$\begin{equation} \sigma_z^{\rm RL} \skew3\tilde{Q}_{\mathrm{R}} \sigma_z^{\rm RL} = -\sigma_z^{\rm RL}+ 2(\mathbbm{1} - Q_+^{\mathrm{T}})\tilde{\cal D}^{-1}. \end{equation} \tag{ C.14 }$$

At  → 0 the inverse operator from the secular matrix has the pole structure, $\tilde {{\cal D}}^{-1} \sim \mathcal {R}/\epsilon$, where the matrix $\mathcal {R}$ is its residue at zero energy.

At this stage it is advantageous to introduce a new bra

$$\begin{equation} \langle \chi_\pm| = \langle \phi_1 |\pm \langle \phi_2|=\langle \phi_1|(\mathbbm{1} \pm Q_+) \end{equation} \tag{ C.15 }$$

and ket

$$\begin{equation} |\chi_\pm \rangle = (\mathbbm{1} \pm Q_+^{\mathrm{T}}) |\phi_1\rangle \end{equation}\noindent \tag{ C.16 }$$

basis in the null spaces Ω^L₀ and Ω^R₀, respectively 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,

$$\begin{equation} \mathbbm{1}=|\eta_\sigma\rangle \langle \chi_\sigma| =|\chi_\sigma\rangle \langle \eta_\sigma|,\end{equation}\noindent \tag{ C.17 }$$

where σ = ±  is summed over.

Now let us look at the singular part of the matrix $\skew3\tilde{Q}_{\mathrm {R}}$,

$$\begin{equation} \sigma_z^{\mathrm{RL}} \skew3\tilde{Q}_{\mathrm{R}}^{\rm sing} \sigma_z^{\mathrm{RL}} = \frac{2}{\epsilon}\mathcal{R} (\mathbbm{1} - Q_+)= \frac{2}{\epsilon} (\mathbbm{1} - Q_+^{\mathrm{T}})\mathcal{R}.\end{equation}\noindent \tag{ C.18 }$$

We note that the matrix $P_-={1\over 2}(\mathbbm{1} - Q_+)$ acts as the projector in the space Ω^L₀, since

$$\begin{equation} \langle \chi_+| P_- = 0\quad \mbox{and} \quad \langle \chi_-| P_-~=~\langle \chi_-|. \end{equation}\noindent \tag{ C.19 }$$

Similar relations hold in the ket space Ω^R₀:

$$\begin{equation} P_-^{\mathrm{T}} |\chi_+\rangle =0 \quad \mbox{and} \quad P_-^{\mathrm{T}} |\chi_-\rangle = |\chi_-\rangle. \end{equation}\noindent \tag{ C.20 }$$

Equivalently, we may write

$$\begin{equation} P_- = |\eta_-\rangle\langle \chi_-|,\,P_-^{\mathrm{T}} = |\chi_-\rangle\langle\eta_-|. \end{equation} \tag{ C.21 }$$

Using these properties, it follows that $\sigma _z^{\mathrm {RL}} \skew3\tilde{Q}_{\mathrm {R}}^{\mathrm { sing}} \sigma _z^{\mathrm {RL}}=\frac {4}{\epsilon } \mathcal {R} P_-=\frac {4}{\epsilon }P_-^{\mathrm {T}} \mathcal {R}$. Multiplying these relations by P₋ from the right, and using the projector property P²₋ = P₋, we obtain the relation

$$\begin{equation} \sigma_z^{\mathrm{RL}} \skew3\tilde{Q}_{\mathrm{R}}^{\rm sing} \sigma_z^{\mathrm{RL}}={4\over \epsilon} P_-^{\mathrm{T}} \mathcal{R} P_-={4\over \epsilon} |\chi_-\rangle\langle\eta_-|\mathcal{R} |\eta_-\rangle\langle \chi_-| \end{equation} \tag{ C.22 }$$

or

$$\begin{equation} \skew3\tilde{Q}_{\mathrm{R}}^{\rm sing} ={4 \mathcal{R}_{--} \over \epsilon} \sigma_z^{\mathrm{RL}} |\chi_-\rangle \langle \chi_-|\sigma_z^{\mathrm{RL}}, \end{equation} \tag{ C.23 }$$

where we defined $\mathcal {R}_{--}\equiv \langle \eta _-|\mathcal {R}|\eta _-\rangle$. Taking into account our definition of the quasiclassical Green's function, equations (29) and (36), we finally obtain that the singular part of the propagator around zero energy at x = x' = x_R takes the form

$$\begin{equation} G(x_{\mathrm{R}},x_{\mathrm{R}};\epsilon) \sim \left(\frac {2 \mathcal{R}_{--}}{\mathrm{i}\epsilon}\right) \hat v^{-1/2}\sigma_z^{\mathrm{RL}}| \chi_- \rangle \langle \chi_-| \sigma_z^{\mathrm{RL}} \hat v^{-1/2}, \end{equation} \tag{ C.24 }$$

where $\hat v$ is the diagonal velocity matrix.

This expression can be compared with the spectral decomposition of the Green function,

$$\begin{eqnarray} &&\fl G(x,x',\epsilon) ={|\psi_0(x)\rangle \langle \psi_0(x')|\over \epsilon}+\sum_{j} \frac{|\psi_{E_j}(x)\rangle\langle \psi_{E_j}(x')|}{\epsilon-E_j} \nonumber\\ &&+ \left( \int_{E_{\mathrm{g}}}^{+\infty} \mathrm{d}E + \int_{-\infty}^{-E_{\mathrm{g}}}\mathrm{d}E \right) \frac{|\psi_{E}(x)\rangle\langle \psi_{E}(x')|}{\epsilon-E}, \end{eqnarray} \tag{ C.25 }$$

where |ψ_E(x)〉 are normalized eigenfunctions of the BdG Hamiltonian, 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 the first term is present if the system contains a Majorana state. One thus concludes that

$$\begin{equation} \left(\frac{2\mathcal{R}^{--}}{\mathrm{i}\hat v}\right)^{-1/2}\sigma_z^{\mathrm{RL}}| \chi_- \rangle = |\psi_{0}(x_{\mathrm{R}})\rangle \end{equation} \tag{ C.26 }$$

is the amplitude of the Majorana particle at the point x_R. The amplitude of the Majorana state at any other point x can be obtained from |ψ₀(x_R)〉 by applying the transfer matrix M(x,x_R).

To conclude, we have shown that if the Majorana number of our system is topologically non-trivial, i.e. ${\cal M}=-1$, then the Green function has the pole at zero energy. Up to a prefactor, the residue at zero energy is then the projector on the one-dimensional linear subspace spanned by the Majorana particle.

In this appendix, we show that the boundary conditions for the Eilenberger matrix Q of a spin–orbit wire terminated at the left end by a (infinitely high) tunnel barrier are equivalent to algebraic relations (60) with the effective matrix Q₋ given by equation (65).

Let us work in the original particle–hole basis where the spinor Ψ takes the form specified by equation (15). We start by rewriting the left boundary conditions (64) in a matrix form. If one defines the reflection matrix

$$\begin{equation} \hat r = \left( \begin{array}{@{}cc@{}} 0 & \mathrm{e}^{\mathrm{i}\phi} \sigma_x^{\mathrm{RL}} \\ {\rm e}^{-{\rm i}\phi} \sigma_x^{\mathrm{RL}} & 0 \end{array} \right) \end{equation} \tag{ D.1 }$$

in the particle space, then $\hat r^*$ is the reflection matrix for holes. Consequently, the spinor Ψ(x_L) satisfies the condition

$$\begin{equation} (\mathbbm{1}-\hat R) \Psi(x_{\mathrm{L}}) = 0, \quad \hat R = \left( \begin{array}{@{}cc@{}} \hat r & 0 \\ 0 & \hat r^* \end{array} \right). \end{equation} \tag{ D.2 }$$

The full reflection matrix $\hat R$ is unitary obeying the particle–hole symmetry, σ^ph_xR^Tσ^ph_x = R. Hence the boundary condition for the bar spinor $\bar \Psi = (\sigma _x^{\mathrm { ph}} \Psi )^{\mathrm {T}}$ takes a similar form

$$\begin{equation} \bar \Psi(x_{\mathrm{L}}) (\mathbbm{1}-\hat R) = 0. \end{equation} \tag{ D.3 }$$

According to the definition (29), the Green function g_e(x,x') inherits the boundary condition right to the tunnel junction from those of the direct product of two spinors, $\Psi (x) \otimes (\bar {\Psi }(x')\, \sigma _z^{\mathrm {RL}} )$. We thus conclude that

$$\begin{eqnarray} &&\bigl(\mathbbm{1}- \hat R\bigr) g(x_{\mathrm{L}},x;\epsilon) = 0,\quad g(x,x_{\mathrm{L}};\epsilon) \bigl(\mathbbm{1} + \hat R\bigr) = 0. \end{eqnarray} \tag{ D.4 }$$

When deriving the second relation, we have taken into account that $\sigma _z^{\mathrm {RL}}\, \hat r\, \sigma _z^{\mathrm {RL}} = -\hat r$. The above conditions are valid for any x > x_L. Sending now x → x_L + 0 and using the relations (B.2) written for the matrix Q_L, one obtains

$$\begin{eqnarray} &&(\mathbbm{1}- \hat R)(\mathbbm{1}-Q_{\mathrm{L}})=0, \quad (\mathbbm{1}+Q_{\mathrm{L}})(\mathbbm{1}-\hat R)=0. \end{eqnarray}\noindent \tag{ D.5 }$$

We see that these boundary conditions are equivalent to the algebraic conditions (60) if one identifies $Q_- = -\hat R$. The normalization $\hat R^2 = \mathbbm{1}$ guarantees that Q₋ belongs to the manifold of the quasiclassical Eilenberger functions. Transforming matrix Q₋ into the Majorana representation (21), we finally obtain the result (65). Note that Pf (Q₋ σ^RL_z) = + 1; thus such a Q-matrix describes a topologically trivial insulator.

The scattering phase ϕ is the parameter of the matrix Q₋ that 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.