Majorana modes and complex band structure of quantum wires

We describe Majorana edge states of a semi-infinite wire using the complex band structure approach. In this method the edge state at a given energy is built as a superposition of evanescent waves. It is shown that the superposition can not always satisfy the required boundary condition, thus restricting the existence of edge modes. We discuss purely 1D and 2D systems, focussing in the latter case on the effect of the Rashba mixing term.


I. INTRODUCTION
The realization of Majorana states as specific excitations of semiconductor quantum wires and of other condensed matter systems has attracted much interest recently.  Experimental evidences of these peculiar states in quantum wires and topological insulators have been presented in Refs. [23][24][25][26] and [27,28], respectively. To assess these evidences, theoretical works are presently focussing on a deeper clarification of the physical properties of Majorana states, including analysis of stability limitations and robustness against disorder and other perturbations. [29][30][31][32][33][34][35] As most characteristic properties, a Majorana state is a zero-energy mode, thus degenerate with the system ground state, that is localized on the system edges or interfaces. Additionally, an energy gap with neighboring states protects the Majorana modes against decoherence. In semiconductor wires, the subject of this work, such zero-energy modes can be induced by the combined action of the following three mechanisms: superconductivity, Rashba spin-orbit coupling and Zeeman magnetic field. Superconductivity introduces particle-hole duality, with excitations lying symmetrically at positive and negative energies with respect to the chemical potential. Rashba spin-orbit coupling introduces helicity by associating motion in space with spin. Thirdly, the Zeeman coupling with an applied magnetic field can be viewed as a tunable interaction with which the system can be placed in different regimes.
The edge character of the Majorana modes suggests a link with the physics of surfaces. Indeed, the similarities with the Schockley states of finite crystals 36 have been pointed out in Ref. [22]. We explore in this work another link with surface physics, namely the description of Majorana modes using the complex band structure of the system. The complex band structure of crystals (or periodic systems in general) is a well established theoretical formalism allowing the generalization from bulk propagating states to states that are localized around defects or interfaces of the crystal. 37 The underlying idea of the generalization is the following: while usual propagating states are described by sustained waves with a real wave number k, there are other solutions with complex wave numbers representing evanescent waves whose amplitude decreases in a given direction (increases in the reversed one).
Waves with complex k cannot be physically realized for infinitely long distances from a given point in both directions due to the divergent behavior in one of them; but they can in a finite domain, or even in a semi-infinite one provided the infinite distance is in the direction of decaying amplitude. The latter situation is precisely the one we expect for edge Majorana modes. While for propagating states the analysis is usually done in terms of the band structure {E α (k)}, where α is the band index and k is an arbitrary real wave number, the complex band structure is presented as {k α (E)}, where for any given energy E there is a set of allowed complex wave numbers k α . Of course, if the energy does not allow propagating states all k α will have an imaginary part. In general, however, propagating and evanescent modes may coexist.
In this work we determine numerically the complex band structure of a semiconductor quantum wire with the above mentioned interactions. Previously, 38 we studied the complex band structure in absence of superconductivity and Zeeman field. Our present numerical approach to Majorana physics complements the analytical methods used in Ref. [21]. We suggest robust numerical algorithms to obtain the complex band structure of 1D and 2D quantum wires in general, discussing the conditions for the existence of Majorana modes in the semi-infinite system. In 2D we focus on the role of the Rashba mixing, giving evidence that the Rashba mixing mechanism hinders the coexistence of two Majorana modes. This is in agreement with the findings of numerical diagonalizations for the finite system. 31 Our results correspond to a continuum Hamiltonian of quantum wires and, from this point of view, they also complement approaches based on tight binding models. 9,14 The analysis in terms of the complex band structure gives the precise conditions for the existence of zeroenergy edge modes, as an alternative to the analysis based on propagating bands. In particular, it clarifies why the topological phase transition from absence to presence of a Majorana mode is signaled by the vanish-ing and reopening of the gap for propagating states with k = 0, as discussed in Ref. [7]. Sections II and III are devoted to the 1D and 2D cases, respectively. The first one is meant as a simple yet illustrative case of the more realistic situation including transverse degrees of freedom. Being well understood, 1-22 the 1D system serves as a benchmark of the complex band structure method.

II. THE 1D CASE
A. Hamiltonian and complex band structure Let us first consider a purely 1D model as in Ref. [21], with motion constrained to be on the x axis. The Hamiltonian reads in this case where the Pauli operators for spin and isospin (in electron-hole space) are represented by the σ's and τ 's, respectively. The successive energy contributions to Eq.
The system exhibits very different physics depending on the field orientation. In 1D the Rashba term in Eq.
(1) is defining a preferred spin direction (y) so that, in abscence of any orbital effects, the physics depends on the field orientation relative to this Rashba-induced direction. Since field orientations along x and z are equivalent we will then consider a field lying on the xy plane, with an azimuthal angle φ B (0 ≤ φ B ≤ 2π). The two limiting cases are the parallel (n =x) and transverse (n =ŷ) orientations with respect to the wire.
The band structure of the infinite system is obtained from Schrödinger equation where the wave function variables are the space coordinate x ∈ (−∞, +∞), the spin η σ ∈ {↑, ↓} and the isospin η τ ∈ {⇑, ⇓}. Spin and isospin basis are taken in z orientation which means, for instance, that |Ψ(x, ↑, ⇑)| 2 is the probability of having at x a quasiparticle with spin pointing along +z (↑) of electron type (⇑). Introducing a state wave number k we now assume the following expansion where the quantum numbers are s σ = ±, s τ = ± and the spin/isospin two-component states fulfill The amplitudes ψ (k) sσsτ in Eq. (3) can be viewed as a set of four complex numbers representing the state for a given wave number k. Substituting Eq. (3) in Eq. (1) we easily get the equations for the state complex amplitudes ψ sσ sτ (we drop label k for convenience) where s = −s. Equation (6) is an algebraic (matrix) equation with nondiagonal couplings of ψ sσ sτ to ψ sσ sτ through pairing and to ψ sσ sτ through the Rashba interaction.
It is important to realize that Eq. (6) has solution only for restricted values of the wave number k. An allowed k which is purely real describes a propagating state while a complex k with nonvanishing imaginary part describes an evanescent state. Finding the physically allowed k's at a given energy E, even numerically, is a non trivial task. Notice, for instance, that the matrix in Eq. (6) is non Hermitian for an arbitrary complex k. While many numerical methods can deal with Hermitian matrices, the methods for non Hermitian operators are scarce. An approach based on the determination of the roots of a high-degree polynomial is possible in 1D. 39 We follow, however, an alternative method allowing an extension to higher dimensions.
We have followed the idea suggested in Ref. [38] to find the allowed k's in a robust and numerically stable way. We arbitrarily choose a particular spin and isospin (s, t) and solve the following linear system if (s σ , s τ ) = (s, t) .
We have found that Eq. (7) is well behaved numerically and always has a solution for any wave number k. In order to find the physically acceptable k's we subsequently look for the zeros of the following function obtained from the solution of Eq. (7), Notice that the state amplitudes and wave numbers fulfilling F 1D (k) = 0 from Eqs. (7) and (8) are also solutions of Eq. (6) and, therefore, they are the physical solutions giving the complex band structure of the system. As an illustrative example, Fig. 1 shows a typical complex band structure. Notice that, contrary to common band representations, we choose the energy as horizontal axis since, within our approach, E is given and the wave numbers are inferred. We only display the positive energy solutions, the negative ones corresponding to the specular image with respect to E = 0. As in Ref. [7] we set our energy units using the Rashba interaction as a reference. Specifically, our length and energy units are The propagating states in Fig. 1a agree with those already discussed in Refs. [7,21]. The evanescent states provide, however, a novel view. We distinguish two types of evanescent modes: those shown in Fig. 1a have a purely imaginary wave number, while those with wave numbers having both real and imaginary parts are shown in panel b). For a given energy, wave numbers always come in sets {k, k * , −k, −k * } or, equivalently, as ±|ℜ(k)| ± i|ℑ(k)|, where ℜ and ℑ stand for real and imaginary parts, respectively. 38 That is, wave numbers come in quartets when both real and imaginary parts are non vanishing (panel b) and in pairs when the wave number is either purely real or purely imaginary (panel a). Figure 1 shows that propagating and evanescent bands can be smoothly connected by the extrema points. In panel a) this connection is clearly seen as horizontal parabolas of opposite curvatures, showing how a purely imaginary wave number transforms into a purely real one as the energy increases. In panel b), ℜ(k) of an evanescent state connects with the extremum of the corresponding propagating band, while the imaginary part disappears when entering the propagating band sector. For a given energy, there are always 8 modes, as corresponds to a quadratic kinetic energy and the 4 × 4 matrix in Eq. (6).
Summarizing this section, we stress that the suggested approach is a robust numerical algorithm allowing the determination of the complex band structure for an arbitrary set of parameters. We will focus next on the description of edge states, in particular zero energy ones representing Majorana modes.

B. Edge states
The complex band structure contains the necessary information on all possible eigenstates. In particular, localized states for which the wave function fades away of some defect or interface are, in general, superpositions of evanescent states. Let us consider a semi-infinite 1D system restricted to x ≥ 0. A priori, we expect the abrupt edge at x = 0 to allow the formation of localized states vanishing at x = 0 and decaying towards x → ∞. Such a state has to be a superposition of evanescent waves having ℑ(k) > 0. Let us assume that the set of such wave numbers is {k}, containingÑ modes for a given energy E. The boundary condition an edge state has to fulfill is then where the C k 's are complex numbers and the ψ (k) sσ sτ 's are the state amplitudes defined in the preceding subsection.
Equation (11) is a set of four linear equations, for (s σ s τ ) = (++, +−, −+, −−), andÑ unknown C k 's. We can immediately conclude that ifÑ < 4 no edge state is possible since there are less unknowns than equations in Eq. (11). That is, if there are less than 4 evanescent modes with ℑ(k) > 0 only the trivial solution C k = 0 is possible, i.e., no physical solution is present. On the contrary, ifÑ > 4 it will always be possible to find a set of nonvanishing C ′ k s. IfÑ = 4 the nontrivial solution will be possible if the system matrix is singular (zero determinant). In that case the condition for the existence of an edge state is For instance, in the case of Fig. 1 the requirement of at least 4 evanescent modes restricts the possibility of an edge mode to energies E < 0.1E SO . We have also computed the determinant of Eq. (12) finding that it vanishes only for E = 0, while it grows linearly for E > 0. Therefore, in this example a zero energy edge state is formed. This is a Majorana mode for which Fig. 2 displays the corresponding density profile, obtained as As expected, the densities vanish at x = 0 and decay with damped oscillations for increasing x. This behavior is in good agreement with previous results. 21 Notice also that because of symmetry the densities remain invariant after spin-isospin inversion. The above results can be taken as examples of the method and its capabilities and we now address the question of how to rationalize the existence of topological phases. It was advanced by Oreg et al. 7 that the signature of an emergent topological phase is the closing and reopening of the k = 0 gap in the band structure of propagating states as one Hamiltonian parameter is varied. Increasing the Zeeman energy in a parallel field geometry the critical parameter for the topological transition is such that zero-energy edge (Majorana) modes are present only when ∆ B > ∆ (c) B . In the transverse field geometry, no zero energy modes are ever present. Let us analyze this scenario from the evanescent band structure point of view. Figure 3 shows the evolution of the wave numbers at zero energy as a function of the Zeeman parameter. The field orientation is along x and, in this case, we always find 8 evanescent modes for any ∆ B . This means that  in this case there are propagating modes at zero energy (dark solid lines appearing for ∆ B > 0.25E SO ). The condition of at least four evanescent modes with ℑ(k) > 0 is fulfilled only when ∆ B < 0.25E SO but, even in this region, no edge state is allowed because of the non vanishing determinant, as shown in panel b).

III. THE 2D CASE
We consider now the generalization of the complex band structure approach to include one extra spatial degree of freedom in the transverse direction to the wire. As before, we assume the wire to be infinite along x, but now motion along y is also possible with the restriction −L y /2 ≤ y ≤ L y /2. This corresponds to a hard wall confinement in y within a square well of length L y . We call this a 2D case although in the literature the name quasi-1D is often used to emphasize that there is confinement in the transverse direction while along the longitudinal one motion is free.

A. Hamiltonian
The generalization of Hamiltonian (1) is where V (y) is the square well of length L y mentioned above. Besides y confinement Hamiltonian (16) includes the transverse kinetic energy and a p y -dependent Rashba term that, following a usual convention, we call the Rashba mixing term. The mixing character of this term can be understood from the fact that the p y operator couples different square well eigenstates. The analysis proceeds in a similar way to the preceding 1D case with the important difference that now the state amplitudes of Eq. (3) become functions of y, Accordingly, the algebraic Eq. (6) becomes a second order differential equation due to the dependence on p 2 y and p y Another difference with respect to the 1D case is the dependence on the magnetic field orientation. Now x and z orientations are no longer physically equivalent due to the orbital currents present in the latter. However, in Eq. (18) we have again assumed a magnetic field in the xy plane with an azimuthal angle φ B since, in this geometry, there are no orbital effects of the field.
From the numerical point of view solving Eqs. (18) for an arbitrary k is much more demanding than the 1D case. We have devised a practical algorithm, 38 introducing a uniform grid in y and an arbitrary matching point y m . The differential equation (18) on all grid points but y m is rewritten using finite difference formulas, with an important caveat: crossing from left to right of the matching point is avoided using non centered finite difference formulas. Another condition is that the state amplitude vanishes at the edges of the y domain. At the matching point y m , instead of Eq. (18), we require the conditions where (s, t) are a pair of arbitrarily chosen spin and isospin labels. In Eq. (19) we have used the notation d (L,R) /dy to indicate derivatives at y m using only left (L) or right (R) neighboring grid points. Of course, a physically valid wave number must correspond to a continuous first derivative of ψ st (y) at the matching point. We thus numerically determine k from the zeros of the function As in the 1D case, when F 2D (k) = 0 is fulfilled the solution from Eq. (19) is equivalent to that of the full Eq. (18). The robustness of the suggested numerical algorithms, both for 2D and 1D, is seen in the fact that for any complex k the F functions can be computed, avoiding singular matrices or ill-behaved numerical problems. 37 Finding the physically allowed k's is then accomplished scanning the complex k plane to determine the zeros of F with the desired accuracy. We have checked that the solution is not affected by the arbitrary choice of state components (s, t) and matching point y m . Besides, since Eq. (18) involves only one variable (y), high numerical precission can be easily obtained using large enough numbers of grid points (N y ≈ 500) and of points for the finite difference derivatives (5 to 11 points).

B. 2D edge states
In 2D the boundary condition for an edge state with x ≥ 0 is similar to Eq. (11), with the state amplitudes changed from scalars to y-functions, We immediately notice that now the boundary condition requires an infinite number of evanescent states {k} due to the infinite number of possible values of y. In practice, a truncation to a restricted set ofÑ evanescent modes has to be done and, accordingly, the condition (21) can be imposed onÑ /4 values of y (the number of modes should be a multiple of 4). Instead of using specific values of y in Eq. (21), we have projected Eq. (21) on the set ofÑ evanescent modes, finding the matrix equation where When the numberÑ of evanescent modes is large enough, the zero eigenvalues of matrix M correspond to Majorana edge states. This Hermitian matrix can be diagonalized numerically with standard methods. For the 1D case of the preceding section the diagonalization of matrix M is equivalent to the approach based on the determinant, Eq. (12), since a zero eigenvalue is associated with a vanishing determinant. In 2D matrix M allows a truncation of the infinite number of evanescent modes. The complete determination of the evanescent bands in 2D is a very tedious task due to the large number of states when varying ∆ B . Instead of finding the equivalent to Fig. 3 for the 2D case, we have followed a different approach to determine the existence of Majorana edge states. Since the topological phase transition is signaled by a complete vanishing of the complex wave number, we have scanned the value of F 2D (k = 0) when varying the Zeeman parameter ∆ B of the Hamiltonian. The regions between nodes of F 2D (k = 0) are the different phases, and we have just studied a few selected values of ∆ B in each region. For a given ∆ B , we look for a set of evanescent modes in the complex-k plane. With these modes we then find matrix M k ′ k and compute its lower eigenvalues. As discussed above, each vanishing eigenvalue is associated with a Majorana zero mode.
In absence of Rashba mixing each transverse mode of the quantum wire behaves as an independent 1D system: there is an onset for the appearance of the n-th band Majorana mode where ε n is the n-th eigenvalue of the transverse quan- tum well, and multiple Majorana modes are possible independently. This physics is contained in the results shown as a dashed line in Fig. 5a where, as could be expected, the derivative discontinuity F 2D (k = 0) has nodes precisely when ∆ B = ∆ (c) B,n . We have also obtained that, in this case, matrix M has n different zero eigenvalues when evaluated for a value of ∆ B such that ∆ (c) B,n+1 . In practice, a fast numerical convergence of these lower eigenvalues to ≈ 10 −9 is easily obtained.
Conspicuous modifications to the above scenario are seen when including the Rashba mixing. Notice that in Fig. 5 we are assuming a regime of strong spin-orbit in which E SO ≈ ε n+1 −ε n . We find that for ∆ B < 0.71E SO , below the first node of Fig. 5a solid line, the matrix M has no zero eigenvalue. In the region between the first and second nodes, 0.71E SO < ∆ B < 2.17E SO , we find clear evidence of a zero eigenvalue, as proved in Fig. 5b. Convergence is, however, rather slow and a large number of evanescent states is required. As a practical scheme, we have included the evanescent states in an ordered way, taking the modulus of the complex wave number as ordering parameter. This criterion gives, for instance, that in Fig. 5b for |k| < 7L −1 SO we include a total of 44 evanescent modes.
In the region 2.17E SO < ∆ B < 3.97E SO , between the second and third nodes of Fig. 5a solid line, the lower eigenvalues of matrix M converge to finite values, as shown in Fig. 5c for ∆ B = 2.5E SO . This behavior is due to the Rashba mixing, for, in its absence, we find two zero eigenvalues in this region. Though numerical, the difference between panels b) and c) of Fig. 5 is quite clear and leads to the conclusion that Rashba mixing hinders the coexistence of two Majorana modes, which is in agreement with the results of Ref. [31] for finite systems. It is remarkable that this hindrance persists even in the semi-infinite system. We focus next on the existence of edge states when ∆ B lies between the third and fourth nodes of F 2D (k = 0); that is, when three zero modes are present without Rashba mixing. As shown in Fig. 6a, one Majorana mode again emerges when enough evanescent states are included in matrix M. This result agrees with the physical scenario found in Ref. 9 using a tight-binding approach; i.e., just a single zero mode is protected in non trivial topological phases. For the sake of discussion, we have artificially weighted the Rashba mixing term by a factor R m in Fig. 6. Convergence with the |k| cutoff is much slower with full mixing (R m = 1) than with a reduced one (R m = 0.2), clearly indicating the relevance of this mechanism.
For completeness, Fig. 7 shows the density distributions corresponding to the Majorana mode of Fig. 5b, using a 2D generalization of Eq. (14). Notice that, on the scale of the figure, the density vanishes for x = 0 and all values of y which is indicating the good convergence of the boundary condition with the number of evanescent modes. Of course, the densities also vanish asymptotically for x → ∞ and, as in the 1D case, remain invariant when inverting spin and isospin.

IV. CONCLUSIONS
We have discussed the physics of Majorana modes in semi-infinite 1D and 2D wires using the complex band structure approach. This formalism provides a natural description of the zero-energy edge modes as superpositions of multiple evanescent waves. The boundary condition of a vanishing wave function at the edge cannot always be fulfilled, this limitation defining the parameter regions (phases) where Majorana modes exist.
The phase transition occurs when for zero energy the complex wave number of an evanescent band vanishes. The phase transition is seen as the emergence of a zero eigenvalue of matrix M, defined in Eq. (23). In 1D the analysis with matrix M is equivalent to the determinant of state amplitudes; while in 2D matrix M allows to monitor the convergence of the lower eigenvalues when the number of evanescent modes is increased. In 2D we find evidence that, in case of a strong Rashba coupling, the Rashba mixing effect hinders the coexistence of multiple Majorana modes in the semi-infinite wire. In agreement with tight-binding models, 9 alternating regions having zero and one Majorana modes are found when increasing ∆ B .