Transport through Majorana nanowires attached to normal leads

This paper presents a coupled channel model for transport in 2D semiconductor Majorana nanowires coupled to normal leads. When the nanowire hosts a zero mode, conspicuous signatures on the linear conductance are predicted. An effective model in second quantization allowing a fully analytical solution is used to clarify the physics. We also discuss the nonlinear current response ($dI/dV$)


Introduction
The theoretical proposal that Majorana modes exist in semiconductor devices [1,2,3] and their subsequent detection in InSb wires [4,5,6] has opened up a new subfield of research on nanostructure properties (see Refs. [7,8,9,10] for reviews). It was originally proposed by Majorana that a massless elementary particle, called a Majorana particle, could exist with the peculiar property of being its own antiparticle. Similarly, a Majorana mode of a semiconductor nanowire is a zero-energy state that remains invariant after charge conjugation. These states are quasiparticle excitations localized on the tips of a finite but long enough wire and they are well separated from the rest of the spectrum of eigenstates by an energy gap.
The existence of Majorana modes in a semiconductor wire requires the presence of the following physical ingredients: a) Zeeman coupling between spin and magnetic field, b) Rashba spin-orbit interaction and c) superconductivity [11,12,13]. The latter can be induced by proximity with a superconductor material and it introduces the concept of electron-hole symmetry [14]. The Rashba spin-orbit interaction is a relativistic effect originating in the quantum well asymmetry in the perpendicular direction to the nanostructure plane. In the present context this interaction introduces chirality by connecting the state of motion with spin. The Zeeman coupling in semiconductors like InAs and InSb is quite large even for relatively low magnetic fields due to the large g factors of these materials. It breaks Kramers degeneracy since the system is no longer time reversal invariant. In this paper we call Majorana nanowire (MNW) a semiconductor nanowire with all three physical effects a), b) and c) mentioned above.
The transport properties in the presence of localized zero modes have been investigated for a normal-superconductor interface [15,16,17]. It has been shown that both for resonant tunneling and for transmission through a quantum point contact a conductance quantization at half-integer multiples of 4e 2 /h is obtained when a zero mode is present at the interface. In this case the superconductor is called topological. In this manuscript we address a related although different geometry, the N/MNW/N structure where N refers to normal contacts in which the pairing and Rashba interactions vanish. We show that the existence of a zero mode is characterized by a unitary Andreev reflection, giving a linear conductance of the N/MNW/N structure equal to e 2 /h. We also find that by increasing the Zeeman coupling parallel to the wire a pronounced dip in the linear conductance appears due to the mixing between channels induced by the Rashba interaction. Similarly to the N/superconductor case [17], the differential conductance has a peak at zero bias in presence of the zero mode, that evolves to a dip when the zero mode is absent. It is also worth stressing that the role of disorder of the N/MNW/N system has been studied in Ref. [18].
Our main contribution in this work is the formalism of the coupled channel model (CCM) for the Bogliubov-deGennes (BdG) Hamiltonian. This formalism can be viewed as an alternative to the methods based on matching of plane waves, or on tight-binding chains [19,20]. It is particularly adapted to the description of spatially smooth potentials and it gives insights on the role of the different physical mechanisms by means of the channel-channel couplings. We present numerical solutions of the CCM equations for a representative case of a 2D MNW based on InAs. In support of our physical interpretations, we also present a simplified effective model allowing a fully analytical solution.

The physical system
The N/MNW/N system is modelled as a 2D channel of transverse dimension L y and with a central region of length L with superconducting and Rashba interactions. These interactions vary smoothly in the longitudinal direction taking constant values ∆ 0 and α 0 in the MNW, and zero in the asymptotic regions of the leads. In addition, potential barriers separate the central MNW from the leads. A sketch of the system and of the x-dependent functions is shown in Fig. 1. The Hamiltonian reads with The x-dependence of the pairing ∆(x), of the Rashba coupling α(x) and of the double barrier V db (x) is modeled by smooth Fermi-like functions with a small diffusivity d; for instance The transverse confinement potential V c (y) is taken simply as an infinite square well with zero potential at the bottom. The chemical potential explicitly appearing in the BdG theory is represented by parameter µ in Eq. (1); while σ and τ are, in a usual notation, the vectors of Pauli matrices acting in spin and isospin (or particle-hole) spaces, respectively. A distinctive characteristic of our model is the continuity of the system parameters with respect to the longitudinal coordinate x, as sketched in Fig. 1b. This would allow us to investigate, for instance, the dependence on the diffusivity d of the transition. In this work, however, we will assume rather steep transitions of the system parameters. The coherent quasiparticle transport is described by the BdG equation where E and Ψ are the quasiparticle energy and wave function, respectively. The latter depends on the position in space (x, y) as well as on the spin and isospin variables (η σ , η τ ), where η =↑, ↓ represents a generic discrete variable with only two possible values, Notice, finally, thatn is assumed to lie in the xy plane. An out-of-plane component would require the addition of orbital magnetic effects not considered in this work.

The coupled channel model
We present in this section the description of transport in terms of channel amplitudes or wave functions obeying a set of coupled differential equations. This description can be viewed as an alternative to the matching of bulk solutions often used in the literature. The CCM is well suited to the problem of a spatially continuous Hamiltonian posed in the preceding section. Before discussing the CCM equations, however, we need to consider the asymptotic solutions (x → ±∞) as they are actually defining the channels themselves.

Asymptotic solutions
Since the pairing and Rashba intensities vanish asymptotically, the BdG Hamiltonian greatly simplifies in those regions, In this limit the eigenstates are spinors pointing in the direction ofn andẑ for spin and isospin. Introducing the quantum numbers s σ = ±1 and s τ = ±1 they read, respectively, where ϕ is the azimuthal angle ofn. The spatial dependence of the asymptotic eigenstates is also analytical, a plane wave in x and a square well eigenfunction in y, φ n (y), with n = 1, 2, . . .. Summarizing, a channel is specified by the quantum numbers (ns σ s τ ) and its wave function reads The propagating or evanescent character of each channel is found when determining its wavenumber k ≡ k nsσsτ . From Eq. (4) the asymptotic BdG energy is where ε n = π 2 n 2 2 2m * L 2 y .
Inverting Eq. (9) it is The channel wavenumber k nsσsτ from Eq. (11) is either real or purely imaginary. These two cases clearly correspond to propagating and evanescent channels, respectively. In conclusion, for electrons (s τ = 1) and holes (s τ = −1) the condition for propagating mode of spin s σ in direction ofn and with transverse state n is

The CCM equations
Assume the following expansion of the full wave function, valid not only in the asymptotic leads but for any arbitrary position, where ψ nsσsτ (x) is a 1D function we call the channel amplitude. Obviously, the channel amplitudes in the asymptotic regions are ψ nsσsτ (x) ∝ exp (±ik nsσsτ x), i.e., propagating or evanescent waves to the right or left directions, depending on the sign of the exponent. The equations fulfilled by the channel amplitudes can be obtained substituting the wave function, Eq. (13), in the BdG equation, Eq. (4), and projecting on a specific channel, The sets of transverse wave functions {φ n }, {χ sσ } and {χ sτ } fulfills proper orthonormality relations. After some straightforward algebra, Eq. (14) leads to where we have introduced the usual anticommutator notation, {p x , α(x)} = p x α(x) + α(x)p x and the bar over an index denotes its opposite value. The set of Eqs. (15) is already a first version of our desired CCM equations. There are three types of contributions to Eq. (15): a) the background terms of channel (ns σ s τ ) are given by the first line, b) the second line contains the coupling terms with channels of the same n but with opposite spins σ or isospins τ to that of the background channel, c) finally, the third line shows the coupling with channels of a different n, the same isospin and arbitrary spin. The physical role played by the different Hamiltonian contributions are clearly seen in Eq. (15). As expected, the superconducting pairing ∆(x) couples electron and hole channels. The two Rashba terms have a markedly different effect regarding the n quantum number. The α(x)p x is diagonal in n, while the α(x)p y is mixing channels with different n's with the selection rules imposed by the square well matrix element n|p y |n ′ . The relevance of the field orientation is also appreciated from Eq. (15). For instance, if the field is along y (ϕ = π/2) the mixing of (ns σ s τ ) and (ns σ s τ ) vanishes.
We end this section mentioning a useful transformation of Eq. (15) that eliminates the linear terms in p x of the background problem. Let us define the transformed channel amplitudeψ where we introduced the dimensionless function Substituting Eq. (16) into Eq. (15) we find The set of Eqs. (18) is very similar to (15), with two important differences: a) the background channel terms have a new contribution quadratic in α(x) which is spin-independent, while the contribution linear in p x is effectively eliminated from this channel, b) the position-dependent phase of the transformation given in Eq. (16) appears explicitly in the coupling with (ns σ s τ ) and (n ′s σ s τ ).

The QTBM
We have solved the set of Eqs. (15) using the quantum-transmitting-boundary formulation of the scattering problem. The reader is addressed to Refs. [22,23] for details on the QTBM. Here we just mention for the sake of completeness the basic underlying ideas. Using a 1D grid Eq. (15) can be discretized with finite-difference formulas for the derivatives. In the asymptotic regions of the leads we impose the analytical solutions of the channel amplitudes where a (i) nsσsτ are the usual incident and reflected amplitudes in lead i. In Eq. (19) we have introduced the lead sign s i , equal to +1 and −1 for the left (i = 1) and right (i = 2) leads, respectively, as well as the position of each lead boundary x i . We have also taken into account the reversed direction of propagation for electrons and holes with the s τ sign. Notice that from Eq. (19) the outgoing coefficient b nsσsτ is expressed in terms of the channel amplitude at the lead boundary, b nsσsτ . (20) The QTBM closed system of linear equations is defined as follows: a) for a grid point x such that x 1 ≤ x ≤ x 2 we impose the discretized version of Eq. (15), b) for a grid point having x < x 1 or x > x 2 we impose Eq. (20). The resulting linear system has as many equations as grid points and it depends only on the set of input coefficients {a (i) nsσsτ }. It is highly sparse and can be numerically solved in an efficient way.
The matrix of transmissions from mode ns σ s τ of lead i to mode n ′ s ′ σ s ′ τ of lead i ′ is given by where the subscript oim, standing for only incident mode, refers to the fact that all incident amplitudes vanish except the one explicitly appearing in the denominator of Eq. (21). For use in the next section, we define a reduced matrix of transmission probabilities where we only discriminate lead and particle type,

Transport in the BdG framework
The description of transport through MNW's can be done with the formalism of transport through normal/superconductor/normal structures. We follow, specifically, the formulation by Lambert et al. [21] for mesoscopic superconductors. For our twoterminal structure, labelled as i = 1, 2 for left and right contacts, the current in terminal i reads where E is the BdG quasiparticle energy. In Eq. (23), J α i (E) andĴ α i (E) are, respectively, the in-going and out-going fluxes in lead i of type α.
The essential ingredients we need to specify in order to use Eq. (23) are the quasiparticle energy distributions f α i (E), the number of propagating modes m α i (E), and the matrix of quantum transmissions P αβ ij (E). The latter two are obtained from the CCM, Eqs. (12) and (22), respectively. The quasiparticle distributions are assumed to be given by Fermi functions where the i-th reservoir chemical potential has been defined as µ i = µ + eV i and kT is the thermal energy. With these inputs the fluxes in Eq. (23) read This formalism fulfills two basic physical conditions: a) vanishing of current for zero bias and b) equality of current in both leads. Indeed, for zero bias all distributions are identical f α i (E) ≡ f (E) and then the sum rule on quantum transmissions, ensures that in-going and out-going fluxes exactly cancel each other. The second condition, I 1 + I 2 = 0, is more subtle; following Lambert [21] we interpret that it actually determines the MNW chemical potential µ, relative to µ 1 and µ 2 . Notice that the potential bias between the two leads is V = V 1 − V 2 and that the MNW chemical potential lies somewhere in the range between the two reservoir chemical potentials, The following practical approach to the BdG transport problem is then suggested: 1) given µ 1 and µ 2 , assume µ = (µ 1 + µ 2 )/2 and solve the BdG-CCM equations for the set {m α i , P αβ ij }; 2) compute I 1 + I 2 ; 3) vary the value of µ and recompute {m α i , P αβ ij } until I 1 + I 2 = 0 is fulfilled. Solving this selfconsistency loop might be a difficult task, however it is not needed when the problem is symmetric with respect to x inversion around the center of the MNW. In this case µ = (µ 1 + µ 2 )/2 is already the solution giving I 1 + I 2 = 0 since the bias V has to be shared symmetrically, V i = s i V /2, where s 1 = 1 and s 2 = −1. Here we shall focus on the symmetric problem, leaving for a future work the analysis on the non symmetric case.

Differential and linear conductances
The differential conductance, defined generically as dI/dV , is one of the most relevant transport properties usually measured in experiments. At zero temperature, the above formalism yields a very simple expression of this quantity because, in this limit, the derivatives of the quasiparticle distribution functions with respect to the bias become Dirac deltas. Of course, this is true only in the symmetric case, when V = 2s i V i .
and, as discussed above, it is dI 2 /dV = −dI 1 /dV . The expression of the linear conductance G can be obtained simply setting the bias to zero in Eq. (29). Using, in addition, the particle hole symmetry we find Equations (29) and (31) are the basic relations of this work. Notice that they contain two qualitatively different contributions to the conductance, a normal transmission, T 0 ≡ (P ++ 12 + P −− 12 )/2, whereby quasiparticle type is conserved; and an Andreev reflection, R A ≡ (P +− 11 + P −+ 11 )/2, with quasiparticle change. Anticipating a result to be discussed below, notice that Eqs. (29) and (31) predict a remarkable phenomenon, a nonvanishing conductance in absence of transmission (T 0 = 0) due solely to Andreev reflection. This occurs when the Majorana nanowire has a zero mode. In this case Andreev reflection is maximal for zero bias, while increasing the bias there is a reduction of R A , i.e., a zero-bias anomaly appears in dI 1 /dV due to the zero mode.

Physical and scaled values of the parameters
The relative strengths of spin-orbit, pairing and Zeeman terms for a given transverse dimension L y are determined by the following scaled dimensionless ratios (scaling is indicated with an s superscript) Notice that for a given set of physical values of α 0 , ∆ 0 and ∆ B different values of the transverse dimension L y will actually correspond to different relative strengths through Eqs. (32-34). Increasing L y the scenario clearly evolves from weak to strong couplings. More specifically, we consider below physical parameters that could represent an InAs-based nanowire [24], m * = 0.033m e , α 0 = 30 meVnm and a pairing gap of ∆ 0 = 0.3 meV. We assume a realistic value of the wire transverse dimension, L y = 150 nm, for which the relative strengths of Rashba and pairing are then α . We also choose the chemical potential, defining our reference energy of the MNW, as µ = 0. Overall, we stress that the complete parameter set is representative of a typical experiment with an InAs-based 2D semiconductor wire.

Linear conductance results
In Figs. 2 to 4 we display the linear conductance calculated from Eq. (31) as a function of the scaled Zeeman value ∆ (s) B for the set of parameters mentioned in the preceding subsection. Figures 2 and 3 correspond to magnetic field in parallel direction to the wire (x) while Fig. 4 is for transverse orientation (y). In Figure 2 we neglected the contribution of the Rashba mixing, i.e., of the terms containing α(x)p y in Eq. (15). Notice that in this situation the linear conductance displays an almost perfect quantization in e 2 /h steps. Increasing ∆ (s) B , small deviations in the form of very narrow spikes can be seen at the beginning of the second and third plateaus. In the first two steps all the conductance is due to Andreev reflection since the normal transmission is negligible. Perfect Andreev reflection is a signal of the existence of a zero mode of the closed system, i.e., a Majorana fermion bound at the interface between the MNW and the normal contacts. This zero mode yields a perfectly quantized conductance in absence of transmission due solely to Andreev reflection, i.e., G = (e 2 /h)R A . The decrease of Andreev reflection for ∆ (s) B > 50 in Fig. 2 can be attributed to the finite size effect that removes the Majorana modes from perfect zero energy, in agreement with the analysis of Ref. [26] for the closed system. This decrease in R A is accompanied by an increase in T 0 , keeping the value of G close to integer multiples of e 2 /h, except at the transition between steps. Figure 3 displays the linear conductance for the same system of Fig. 2, but now including the full Rashba interaction. A conspicuous difference with Fig. 2 is that the conductance deviates from the simple staircase behaviour, with a broad conductance dip appearing at the end of the first plateau. This dip is due to a magnetic instability precluding the formation of two zero modes due to a repulsion between modes induced by Rashba mixing [26]. The effect of this mechanism on the linear conductance is remarkable, with the prediction of a reduced conductance due to a large reduction of Andreev reflection. This anomalous behaviour of the conductance at the end of the conductance plateau also appears in higher plateaus, as seen in Fig. 3 for the second and third plateaus. Notice, however, that the finite size effect mentioned above transforms the conductance dips of the higher plateaus in a strongly oscillating behaviour. We have checked that the formation of the conductance dips due to the instability of multiple zero modes is even more robust with higher values of the pairing gap and Rashba strengths, and that it is also robust against variations of the barriers between the normal contacts and the MNW. Figure 4 shows the evolution of the linear conductance with ∆ (s) B for field along the transverse direction y. For this orientation of the field the physics changes completely, since now it is the Andreev reflection that vanishes and the conductance is due to the normal transmission. Only small peaks in R A can be seen in the transition between plateaus. The vanishing of R A is due to the absence of zero modes of the closed MNW for magnetic fields along y [26]. A similar orientation anisotropy has been seen in experiments with cylindric InSb nanowires [4,5]. No conductance dips are observed in Fig. 4 but there are many spikes due to resonant transmission through the doublebarrier potential V db (x). The separation between spikes is very small due to the dense distribution of quasibound states for such a long system L = 3 µm. From this point of view, it is still more remarkable the fact that for field along x the presence of a zero mode washes the spike oscillations and yields a consistent maximal conductance in some regimes.

Nonlinear conductance
The nonlinear conductance obtained with Eq. (29) is shown in Fig. 5 as a function of the applied bias. We have taken some selected values of ∆ (s) B from Fig. 3, corresponding to vanishing bias, and explored the variation with V . As in the preceding subsection we define a scaled bias taking the transverse confinement as reference, i.e., V (s) = (em * L 2 y / 2 )V . For ∆ (s) B = 9 there is a narrow peak in dI 1 /dV at zero bias. This zero bias anomaly is reflecting the existence of a zero mode in the MNW. Increasing the Zeeman coupling the peak broadens, becoming a flat distribution. For ∆ (s) B = 19, corresponding to the conductance dip of Fig. 3, the zero bias peak changes to a zero bias minimum. The existence of a zero bias anomaly in the presence of zero modes has been discussed before in systems with a superconductor contact, experimentally in Refs. [4,5] and theoretically in Refs. [15,16,17]. Our results prove that a similar behaviour is to be expected in N/MNW/N structures.

Density distributions
The density distributions, defined as |ψ nsσsτ (x)| 2 , are shown in Fig. 6 for two values of ∆  Fig. 3. As expected, the upper panel shows that the incident unitary density couples with an edge mode of the MNW. The density profile localized at the edge and decaying towards the interior has exactly the same shape found in calculations of zero modes of closed MNW's [26]. Perfect Andreev reflection in this situation consists in the total reflection in the conjugate channel and, therefore, no quantum interference is observed in the left contact. As mentioned before, this occurs due to the presence of the zero mode in the MNW and it allows unit conductance without any transmission at all between left and right contacts.
The lower panel of Fig. 6 shows a qualitatively different behaviour. The beating pattern in the left contact is indicating that full reflection occurs now in the same channel of incidence, with a strong interference between incident and reflected waves The density at the edge of the MNW is more irregular and extends farther towards the interior than in the upper panel. The physical interpretation is clear: for this Zeeman intensity the edge mode of the MNW lies at a nonzero energy, this causing normal reflection, as opposed to the Andreev reflection of the upper panel.
To complete the calculation we need to determine the MNW Green's functions. They read [16] Here, and the self-energy matrices are given by Substituting Eqs. (40)-(45) into Eq. (39) and using current conservation we obtain From Eq. (46) we note that at T = 0 the linear conductance finally reads For zero energy Majoranas ε M = 0 and then Eq. (46) yields G = e 2 /h. This result nicely agrees with our interpretation of Fig. 3 attributing maximal conductance to the zero mode.

Conclusions
We have presented the formalism of transport in a N/MNW/N structure based on the coupled channel model. This formalism yields a transparent interpretation of the coupling between channels induced by the relevant physical mechanisms of the problem. Namely, the confinement, Zeeman, Rashba and superconducting interactions. We have considered a 2D structure and in-plane magnetic fields, although the formalism can be extended to consider more spatial dimensions and different geometries.
The coupled-channel-model equations have been solved using the quantumtransmitting-boundary algorithm for a set of parameters representative of an InAs nanowire. The existence of a zero mode in the MNW is characterized by a perfect Andreev reflection, whereby an incident channel is totally reflected in its antiparticle conjugate one. For a single zero mode the linear conductance takes the maximal value e 2 /h due solely to Andreev reflection, without any quantum transmission from left to right contacts. For increasing values of the Zeeman coupling along the wire, a conspicuous dip in the linear conductance is predicted due to repulsion between Majoranas. This repulsion originates in the Rashba mixing between channels. On the contrary, for Zeeman coupling along y the Andreev reflection vanishes, with the possible exception of a small region close to the transition between plateaus. When the zero mode is absent, the linear conductance has narrow spikes as a function of the Zeeman coupling.
The differential conductance signals the presence of the zero mode with a peak at zero bias. The zero bias peak evolves to a dip when the MNW zero mode is absent. Finally, we have also discussed an effective model in second quantization confirming the physical interpretation in terms of Majorana modes. The coupled channel model presented in this work can be used to investigate other scenarios like, e.g., non-symmetric barriers or sequential MNW's. Work along these lines is in progress.