Leakage of Majorana bound states in an inhomogeneous topological nanowire

The search of Majorana fermions [1] have attracted increased attention in recent years in condensed matter physics [2–6], due to their fundamental interest and potential application in topological quantum computing [7–13]. In solid states, the evidence of Majorana bound states (MBSs) has popped up at the ends of semiconductor–superconductor hybrid nanowires [14–16], chains of single atoms [17], iron-based superconductors [18–20], topological insulators [21, 22], and the quantum anomalous Hall insulator [23]. Of all the Majorana platforms, the one that has attracted the most excitement is the system based on semiconductor nanowire with Rashba spin–orbit interaction and in proximity to an s-wave superconductor [2, 15, 24]. Among several smoking gun signals for MBSs, the standout is the so-called zero-bias conductance peak, which has been observed in the charge transport measurements of both nanowires [25, 26] and spin chains [27]. However, the zero-bias signature could also arise from other non-Majorana mechanisms [24, 28]. So far, the existence of Majorana quasiparticles has been one of the most important issues in the study of topological superconductors.

One possibility to check the existence of MBSs is to look for signs that the MBSs are hybridizing with each other. For a finite nanowire, the wave functions of MBSs decay exponentially along the wire and their hybridization produces a finite Majorana energy splitting _M ∼ exp(−L/ξ) [29], where L is the length of the nanowire and ξ is the coherence length of the superconductor. The strength of Majorana energy splitting could be tuned periodically by the Zeeman field or the chemical potential of the nanowire. The exponential suppression of energy splitting with increasing wire length is observed in the InAs nanowire segments with epitaxial aluminium, which forms a proximity-induced superconducting Coulomb island [30]. The oscillation of the Majorana energy splitting as a function of the magnetic field is also demonstrated in the experiment [30]. However, the energy splitting indicates a damping oscillation with increase in magnetic field, which contradicts to the theoretical prediction [29, 31]. Recently, a possible explanation for this observation is proposed by considering a position-dependent spin–orbit coupling [32]. This implies that the inhomogeneity in the nanowire could induce many significant physical consequences [33].

Theoretically, since MBSs were proposed to be realized in a semiconductor nanowire-superconductor hybrid system, numerous studies have analyzed the model Hamiltonian by using the tight-binding model and direct numerical diagonalization [34–37]. Up to now, the solution of the continuum Bogolyubov-de-Gennes (BdG) equation about non-zero energy Majorana mode, is based on the variational calculation or small expansion [31]. In the long wire length limit, the approximation is appropriate. However, it is difficult to judge its accuracy in other complex regimes. The exact solution of the continuum Hamiltonian for MBSs in the finite-size case is expected to provide additional insights into the nature of MBSs. An exact solution for the continuum BdG equation of the superconductor–semiconductor hybrid nanowire system is still missing.

Motivated by this, here we present an exact solution of the continuum BdG equation for MBSs with a finite splitting energy in a topological superconductive nanowire by means of the trial wave function method [38]. The open boundary condition at wire ends is taken in the calculation. This method could give the full energy spectra of the BdG Hamiltonian and the corresponding wave functions. The paper is organized as follows. In section 2 we introduce the BdG model for the MBSs and provide the calculation procedure within an exact treatment. In section 3 we discuss the dependence of the splitting energy on the system parameters in a homogeneous nanowire. We extend the discussion to the inhomogeneous case in section 4. The spatial distribution of MBS wave functions is investigated after applying a local gate voltage. Finally, a summary is presented in section 5.

We focus on one of the most promising system for generating MBSs, which is by inducing superconductivity on a semiconductor nanowire with Rashba spin–orbit coupling. Under a proper application of external Zeeman field, the nanowire resembles a one-dimensional topological superconductor hosting MBSs at its two ends. Following the earlier work, the one-dimensional BdG Hamiltonian for the present system can be written as [14]

$$\begin{equation}{H}_{\mathrm{B}\mathrm{d}\mathrm{G}}=\left[-\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{\partial }_{x}^{2}-\mu \right]{\tau }_{z}+{V}_{\mathrm{z}}{\sigma }_{z}+i\alpha {\partial }_{x}{\sigma }_{y}{\tau }_{z}+{\Delta}{\tau }_{x},\end{equation} \tag{ 1 }$$

where V_z is the effective Zeeman field, σ_i and τ_i (i = x, y, z) are Pauli matrices acting on the spin and particle–hole spaces, μ is the system chemical potential strength, α represents the spin–orbit coupling, and Δ is the proximity-induced superconductor gap in the wire. Here we concentrate on the exact one-dimensional case and study the finite-size effects.

We study the non-zero-energy solution of MBSs by employing the trial wave function method. Due to the finite size of the system, the non-zero-energy solution of MBSs could be obtained, and their spinor wave functions do not possess a standard form as the zero-energy mode with Ψ = (ψ, iσ_yψ*) any more. In order to get the non-zero-energy solution of MBSs, we need to solve the BdG equation H_BdGΨ(x) = EΨ(x) with vanishing boundary conditions, where E is the eigenenergy and Ψ(x) is the corresponding wave function. We solve the BdG equation by using the trial wave function of the following form

$$\begin{equation}{\Psi}\left(x\right)=\sum _{n}{a}_{n}{u}_{n}\left(x\right)=\sum _{n}{a}_{n}{\rho }_{n}\enspace {\mathrm{e}}^{{z}_{n}x},\end{equation} \tag{ 2 }$$

where u_n(x) are the eigenstate of the equation H_BdGu_n(x) = Eu_n(x) and are assumed to have the form ${\rho }_{n}\enspace {\mathrm{e}}^{{z}_{n}x}$, a_n are the corresponding coefficients. The z_n of a non-trivial solution must satisfy the quartic equation that can be obtained from the secular equation about z_n

$$\begin{align}\hfill & {\left[{\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}^{2}+\mu \right)}^{2}+{\alpha }^{2}{z}^{2}+{{\Delta}}^{2}-{V}_{\mathrm{z}}^{2}-{E}^{2}\right]}^{2}\hfill \\ \hfill & \quad +\enspace 4{\alpha }^{2}{z}^{2}\left({E}^{2}-{{\Delta}}^{2}\right)-4{E}^{2}{V}_{\mathrm{z}}^{2}=0.\hfill \end{align} \tag{ 3 }$$

This is a strict relationship between an eigenenergy and its wave function. For a semi-infinite nanowire case, the Majorana mode is exact zero energy and its wave function has been investigated in detail [31]. The condition that Majorana zero modes appear in the nanowire is ${V}_{\mathrm{z}}^{2}-{{\Delta}}^{2}-{\mu }^{2}{ >}0$ [15].

After solving the equation (3), we can obtain all values of z_n for a given E. For each z_n, the corresponding u_n(x) without renormalization is given by

$$\begin{equation}{u}_{n}\left(x\right)={\rho }_{n}\enspace {\mathrm{e}}^{{z}_{n}x}=\left(\begin{pmatrix}\hfill {u}_{n1}\hfill \\ \hfill {u}_{n2}\hfill \\ \hfill {u}_{n3}\hfill \\ \hfill {u}_{n4}\hfill \end{pmatrix}\right){\mathrm{e}}^{{z}_{n}x},\end{equation} \tag{ 4 }$$

and the four components are

$$\begin{align*}\hfill {u}_{n1}& =2\alpha {z}_{n}\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+\mu +{V}_{\mathrm{z}}\right)\hfill \\ \hfill {u}_{n2}& ={\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+\mu \right)}^{2}-{\alpha }^{2}{z}_{n}^{2}-{\left(E-{V}_{\mathrm{z}}\right)}^{2}+{{\Delta}}^{2}\hfill \\ \hfill {u}_{n3}& =-2{C}_{n}\alpha {z}_{n}\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+\mu -{V}_{\mathrm{z}}\right)\hfill \\ \hfill {u}_{n4}& ={C}_{n}\left[{\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+\mu \right)}^{2}-{\alpha }^{2}{z}_{n}^{2}-{\left(E+{V}_{\mathrm{z}}\right)}^{2}+{{\Delta}}^{2}\right]\hfill \end{align*}$$

with

$$\begin{equation}{C}_{n}=\frac{{\left(\frac{{\hslash }^{2}{z}_{n}^{2}}{2{m}^{{\ast}}}+\mu +E\right)}^{2}+{\alpha }^{2}{z}_{n}^{2}-{V}_{\mathrm{z}}^{2}-{{\Delta}}^{2}\left.\right]}{\left[2{\Delta}\left(\frac{{\hslash }^{2}{z}_{n}^{2}}{2{m}^{{\ast}}}+\mu -{V}_{\mathrm{z}}\right)\right]}.\end{equation} \tag{ 5 }$$

One can prove that C_n = 1 for the exact zero-energy mode E = 0.

To determine the energy of in-gap bound states, the boundary conditions of a wave function for a wire with finite length are needed. At the boundary x = ±L/2 (L is the wire length), we impose boundary conditions such that the wave function vanishes at the ends, Ψ(±L/2) = 0. Combining equations (2) with (4), we can get eight equations about a_n (n = 1, 2..., 8) from boundary conditions ${\Psi}\left({\pm}L/2\right)={\sum }_{n=1}^{8}{a}_{n}\enspace {\mathrm{e}}^{{\pm}{z}_{n}L/2}{\rho }_{n}=0$. It is required that the coefficient matrix about a_n has a vanishing determinant and the secular equation about E becomes

$$\begin{equation}\mathrm{D}\mathrm{e}\mathrm{t}\left({\mathbf{M}}^{{\pm}}\right)=\left\vert \begin{matrix}\hfill {u}_{11}{t}_{1}^{{\pm}}\hfill & \hfill {u}_{21}{t}_{2}^{{\pm}}\hfill & \hfill {u}_{31}{t}_{3}^{{\pm}}\hfill & \hfill {u}_{41}{t}_{4}^{{\pm}}\hfill \\ \hfill {u}_{12}\hfill & \hfill {u}_{22}\hfill & \hfill {u}_{32}\hfill & \hfill {u}_{42}\hfill \\ \hfill {u}_{13}{t}_{1}^{{\pm}}\hfill & \hfill {u}_{23}{t}_{2}^{{\pm}}\hfill & \hfill {u}_{33}{t}_{3}^{{\pm}}\hfill & \hfill {u}_{43}{t}_{4}^{{\pm}}\hfill \\ \hfill {u}_{14}\hfill & \hfill {u}_{24}\hfill & \hfill {u}_{34}\hfill & \hfill {u}_{44}\hfill \end{matrix}\right\vert =0,\end{equation} \tag{ 6 }$$

where ${t}_{m}^{+}=\mathrm{tanh}\enspace \frac{{z}_{m}L}{2}$, ${t}_{m}^{-}=c\enspace \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\frac{{z}_{m}L}{2}$ and z_m are four roots with non-negative real part. Details are given in the appendix A. Det(M^±) = 0 corresponds to the particle or hole solution of E, respectively. We determine E is the eigenenergy or not by judging whether at least one of Det(M^±) is zero or not. Once an energy E meets this condition, it is considered as an eigenenergy of the system. The trial wave function method could give all solutions of the BdG Hamiltonian in principle within a nonperturbative treatment.

Because of the complexity of the Hamiltonian, it is difficult to obtain the analytical expression of Majorana energy splitting. The energy spectrum and the wave functions are solved numerically from equation (6). However, for strong enough Zeeman splitting, the explicit formula about Majorana energy splitting is presented by trial function method [as shown in equation (13)]. It is noted that the numerical methods in general refer to the tight-binding calculation. Here the trial wave function method is used to solve the BdG equation with a set of proper boundary conditions, which usually refers to the analytical treatment.

In the following, we present and discuss the numerical results on the energy spectra and the corresponding wave functions of the semiconductor–superconductor hybrid system in the homogeneous and the inhomogeneous case respectively. We derive an analytic formula of the splitting energy of Majorana modes for a finite effective p-wave superconducting nanowire.

In figure 1, we present the splitting energy of Majorana modes as a function of Zeeman field V_z for different wire lengths with the chemical potential μ = 0.2 meV. In addition to the energy levels of MBSs, the full low-energy spectra are also illustrated in figure 1. Two lowest-energy spectra obtained by using tight binding method are plotted for comparison. When performing the tight-binding method (TB), we usually need to firstly discretize the continuum model into a lattice model, where only the hopping between the nearest neighbor sites is considered. Therefore, the energy spectrum and eigenstates obtained from the tight-binding method are the approximated results. Differently, no approximation is considered when performing the trial-wave functions method. In this sense, the trial wave function method is better than TB. When V_z exceeds a critical value ${V}_{\mathrm{z}}^{\mathrm{c}}$ (${V}_{\mathrm{z}}^{\mathrm{c}}=\sqrt{{{\Delta}}^{2}+{\mu }^{2}}$), the energy band experiences a closing and reopening process, corresponding to the topological phase transition. For ${V}_{\mathrm{z}}{ >}{V}_{\mathrm{z}}^{\mathrm{c}}$, a set of Majorana energy mode appears near the zero energy. It is shown that the energy splitting of Majorana modes is oscillatory and the splitting strength becomes stronger with the increase of Zeeman field. For a longer nanowire, the finite-size effect becomes weaker and the splitting energy is suppressed, while the oscillation period becomes smaller. Importantly, the oscillation of zero mode splitting is always accompanied by the oscillation of other in-gap states. Different from the MBSs, they are prohibited by the topology. These in-gap states and their oscillatory behaviors are a manifestation of the finite-size effect. Our numerical results indicate that the oscillation magnitude of band gap is comparable to that of zero mode splitting for relatively large Zeeman field.

Zoom In Zoom Out Reset image size

In figure 2, we demonstrate the splitting energy of Majorana modes as a function of the chemical potential for different wire lengths with applied external Zeeman field V_z = 1.0 meV. Two lowest-energy spectra obtained by using tight binding method are also plotted. The energy gap closes at two critical points ±μ_c (${\mu }_{\mathrm{c}}=\sqrt{{V}_{\mathrm{z}}^{2}-{{\Delta}}^{2}}$), and the MBSs appear in the topological regime (−μ_c, μ_c). Similar to the case of varying Zeeman field, the splitting energy of Majorana mode also shows an oscillation behavior and the oscillatory magnitude increases with the increase of the chemical potential. For μ < −μ_c, no discrete energy levels are found in the energy gap. However, in another topological trivial phase with μ > μ_c, there are several discrete energy levels appear in the gap, as shown in figure 2. For a relatively short wire, even the zero-energy mode appears and it also oscillates as a function of chemical potential. The zero-energy modes in the trivial phase are induced by the finite-size effect. When we increase the chemical potential, these modes have a different behavior to that of MBSs in the topological phase. They are not topologically protected and easier to be destroyed by disorder. For the case of longer wire, the oscillation of zero energy mode in the topological phase becomes weaker and its period is smaller. In the topological trivial phase, the zero energy mode disappears gradually, and the discrete energy modes are away from zero energy further with the increase of the chemical potential of the nanowire.

Zoom In Zoom Out Reset image size

In the strong Zeeman field regime that V_z ≫ m*α²/ℏ², Δ, the spins nearly polarize and one can project the original Hamiltonian into a simpler one-band problem. To leading order, one obtains an effective p-wave like Hamiltonian [9, 38–40]

$$\begin{equation}{H}_{\text{eff}}=\left(\begin{pmatrix}\hfill -\frac{{\hslash }^{2}}{2{m}^{{\ast}}}\frac{{\partial }^{2}}{\partial {x}^{2}}-{\mu }_{\text{eff}}\hfill & \hfill -i\hslash {{\Delta}}_{\text{eff}}\frac{\partial }{\partial x}\hfill \\ \hfill -i\hslash {{\Delta}}_{\text{eff}}\frac{\partial }{\partial x}\hfill & \hfill \frac{{\hslash }^{2}}{2{m}^{{\ast}}}\frac{{\partial }^{2}}{\partial {x}^{2}}+{\mu }_{\text{eff}}\hfill \end{pmatrix}\right),\end{equation} \tag{ 7 }$$

where μ_eff = μ + V_z/2 and the effective p-wave pairing field is Δ_eff = αΔ/V_z. Assuming the trial wave function

$$\begin{equation}{\psi }_{n}\left(x\right)=\left(\begin{pmatrix}\hfill {a}_{n}\hfill \\ \hfill {b}_{n}\hfill \end{pmatrix}\right){\mathrm{e}}^{{z}_{n}x}\end{equation} \tag{ 8 }$$

and substituting it into the eigenequation, the secular equation about z_n is given by

$$\begin{equation}{E}^{2}-{\left(\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+{\mu }_{\text{eff}}\right)}^{2}+{\hslash }^{2}{{\Delta}}_{\text{eff}}^{2}{z}_{n}^{2}=0.\end{equation} \tag{ 9 }$$

The wave function Ψ(x) of the eigenstates can be considered as a linear combination of ψ_n(x),

$$\begin{equation}{\Psi}\left(x\right)=\sum _{n}{c}_{n}\left(\begin{pmatrix}\hfill {a}_{n}\hfill \\ \hfill {b}_{n}\hfill \end{pmatrix}\right){\mathrm{e}}^{{z}_{n}x},\end{equation} \tag{ 10 }$$

where

$$\begin{equation}\frac{{a}_{n}}{{b}_{n}}=-\frac{\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{n}^{2}+{\mu }_{\text{eff}}-E}{i\hslash {{\Delta}}_{\text{eff}}{z}_{n}}.\end{equation} \tag{ 11 }$$

By using the vanishing boundary condition ψ(x = ±L/2) = 0, we finally arrive at the transcendental equation for the bound state energies

$$\begin{equation}\frac{\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{1}^{2}+{\mu }_{\text{eff}}-E}{{z}_{1}\enspace \mathrm{tanh}\left(\frac{{z}_{1}L}{2}\right)}=\frac{\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{z}_{2}^{2}+{\mu }_{\text{eff}}-E}{{z}_{2}\enspace \mathrm{tanh}\left(\frac{{z}_{2}L}{2}\right)},\end{equation} \tag{ 12 }$$

where z₁ and z₂ are the roots with positive real part. After some algebra, we obtain an analytical expression of the Majorana energy splitting

$$\begin{equation}{E}_{\mathrm{M}}\approx {\mu }_{\text{eff}}\sqrt{\frac{2{\mu }_{\text{eff}}}{{m}^{{\ast}}{{\Delta}}_{\text{eff}}^{2}}}\enspace {\mathrm{e}}^{-\frac{{m}^{{\ast}}{{\Delta}}_{\text{eff}}L}{\hslash }}\enspace \mathrm{sin}\left(\frac{\sqrt{2{\mu }_{\text{eff}}{m}^{{\ast}}}}{\hslash }L\right).\end{equation} \tag{ 13 }$$

The appearance of the sine function with a function of L, V_z and μ as its arguments can account for the oscillation behaviors of the Majorana energy splitting. The factor in front of the sine function, which is an increasing function of μ or V_z, gives rise to an increasing oscillation amplitude. One can also notice an exponential damping of the oscillation amplitude as a function of L.

Once an eigenenergy of the system is found, its corresponding wave function and probability density distribution in real space could be obtained. Since Hamiltonian equation (1) is written down in terms of the Nambu spinor $\left({c}_{{\uparrow}}\left(x\right),{c}_{{\downarrow}}\left(x\right),{c}_{{\downarrow}}^{{\dagger}}\left(x\right),-{c}_{{\uparrow}}^{{\dagger}}\left(x\right)\right)$, the wave function of an eigenstate takes the form Ψ(x) = (ψ_↑(x), ψ_↓(x), ϕ_↓(x), ϕ_↑(x)). In figure 3, we present the probability density |Ψ_σ(x)|² = |ψ_σ(x)|² + |ϕ_σ(x)|² (σ = ↑, ↓) and |Ψ(x)|² = |Ψ_↓(x)|² + |Ψ_↑(x)|² of the wave function along the wire for different Zeeman fields with a fixed chemical potential μ = 0.2 meV and different chemical potentials with a fixed Zeeman field V_z = 0.8 meV. The nanowire we considered has a length L = 2 μm. From the above discussions, the finite-size effect is supposed to be more weaker. Near the transition point ${V}_{\mathrm{z}}^{\mathrm{c}}$, the wave function for MBSs is mainly localized at both ends of the wire, and there is no oscillatory overlap between two well-separated wave packets. In the case where V_z is small, two spin parts |Ψ_↓(x)|² and |Ψ_↑(x)|² contribute comparably to the |Ψ(x)|², as shown in figure 3(a). With the increase of V_z, the main wave packets move along the wire gradually and the distribution indicates an oscillatory behavior. For a strong enough V_z, the distribution of the MBSs is not restricted to the wire ends any more and |Ψ(x)|² mainly comes from |Ψ_↓(x)|² [see figures 3(c) and (d)]. Majorana modes could nearly polarize in a specified direction [9], here, the z direction. Other Majorana spin related properties have been investigated [41]. Although a Majorana mode is stable in the topological region, the wave function of a Majorana mode is almost equally distributed in the wire for large V_z, as indicated in figure 3(d). In this case, the topologically protected Majorana mode is not bound states any more. In the figures 3(e)–(h), we vary the chemical potential of the nanowire. The applied Zeeman field V_z = 0.8 meV is large enough to show that the spin component of wave density |Ψ_↑(x)|² is suppressed. At the lower critical point μ = −μ_c, the wave function of MBSs is localized at the ends of a nanowire and no oscillatory overlap is formed. As we tune the chemical potential away from this critical point and into the topological regime, the overlap of wave functions of two end MBSs becomes oscillatory. Near two critical points μ = ±μ_c, the Majorana zero mode shows quite different behavior. For μ ∼ μ_c, there are two imaginary roots and the Majoarana wave function displays the envelope oscillation. While for μ ∼ −μ_c, all the roots are real and the wave function exponentially decays along the wire.

Zoom In Zoom Out Reset image size

The present trial wave function method also allows us to extend the discussion to the inhomogeneous nanowire case, which is expected to provide new evidence for the Majorana existence. In this section we consider an inhomogeneous case where a nanowire has two segments with different chemical potentials. The proposed scheme is shown in figure 4(a). The lengths of two segments are L₁ = 1.5 μm, L₂ = 0.5 μm respectively, and are assumed to be much larger than any decay length. This assumption means the Majorana energy splitting is much smaller than the induced superconducting gap. To solve the BdG equation of the entire nanowire, we invoke continuity of the wave functions and their first derivative, i.e. Ψ(x = 0⁺) = Ψ(x = 0⁻), ∂_xΨ(x = 0⁺) = ∂_xΨ(x = 0⁻), and impose vanishing boundary conditions at two ends, Ψ(x = −L₁) = Ψ(x = L₂) = 0. The eigenenergies and their wave functions of the inhomogeneous nanowire can be found in various regimes by matching the boundary conditions.

Zoom In Zoom Out Reset image size

In figure 5, the low-energy spectra of full system are plotted as a function of the chemical potential μ₂ controlled by a local gate voltage applied on the shorter segment. The longer segment is in topological superconducting phase. In the topological nontrivial regime of the shorter segment, it is shown in figure 5 that the MBS energy is periodically modulated by the local gate voltage. The barrier region is formed in the shorter segment of the nanowire for −μ₂ > −μ₁, while a well region is formed in the shorter segment in the case of for −μ₂ < −μ₁. For the potential well case with the chemical potential exceeding the upper topological transition point ${\mu }_{2}^{\mathrm{c}}$ (${\mu }_{2}^{\mathrm{c}}=\sqrt{{V}_{\mathrm{z}}^{2}-{{\Delta}}^{2}}$), there still exist some oscillatory near-zero-energy modes, while such near-zero-energy states are absent in the barrier case where the chemical potential of the shorter segment is smaller than the lower topological transition value $-{\mu }_{2}^{\mathrm{c}}$. Figure 5 suggests that the Majorana energy oscillation can be tuned by locally modulating the chemical potential in one part of the nanowire.

Zoom In Zoom Out Reset image size

We present the spatial probability density distribution of the Majorana wave function by tuning the local gate voltage in figure 6. It is illustrated that the Majorana wave function is almost spatially localized in the topological superconductor region for the barrier case [see figures 6(a) and (b)]. With the increase of the local gate voltage, the entire nanowire is tuned into the topological regime and the MBSs are localized at two wire ends asymmetrically due to the chemical potential difference of two segments [see figure 6(c)]. As the local gate voltage increases further and the shorter segment is driven out of the topological regime, the Majorana state could also penetrate into the well region, as shown in figures 6(d)–(f). There exists a near-zero-energy state in the well region and it could merge with the Majorana state. The wave functions of the lowest-energy bound states in these two regimes indicate different behaviors, but they could coexist in the nanowire. As shown in right panels of figure 6, it is interesting that both spin components of the trivial Andreev bound states in the potential well region account for a significant proportion, while MBSs in the topological region are almost fully spin-polarized.

Zoom In Zoom Out Reset image size

We now proceed to investigate the dependence of the leakage of MBSs on the local gate voltage μ₂ in the gated region. In figure 7, we plot total probability P₂ of the lowest eigenstate in the gated region of the nanowire as a function of the gate voltage μ₂. We see that the leakage of the MBSs depends significantly on the gate voltage. For a low enough gate voltage case that ${\mu }_{2}{< }-{\mu }_{2}^{\mathrm{c}}$, the wave functions of the MBSs are almost localized in the topological region. In this case, the probability densities |Ψ|² in the gated region are exponentially suppressed with the decrease of μ₂. For the case that the gated region is in the topological regime, i.e., $-{\mu }_{2}^{\mathrm{c}}{< }{\mu }_{2}{< }{\mu }_{2}^{\mathrm{c}}$, the MBSs are well located in the two wire ends and P₂ equals to 1/2 approximately. As μ₂ increases further and the gated region is in the trivial phase again, the probability leakage magnitude indicates strong oscillation and does not decay even for the strong enough well, as shown in figure 7. Indeed, for ${\mu }_{2}{ >}{\mu }_{2}^{\mathrm{c}}$ the near-zero-energy Andreev bound states in the trivial region are expected to live, as proved in figure 2 and hybridize with MBSs in the topological region. When an inhomogeneous topological nanowire with part of it being topological trivial, strictly speaking, we can not tell whether a lowest-energy state is topological or not. Physically, we can treat it as the coalescence of a topological state and an Andreev bound state. However, the appearance of these Andreev bound states in topological trivial region could strongly influence the local measurement of Majorana properties [24]. In the experiments, the tunnel probe is attached to the gated part of the nanowire and is thus predominantly coupled to the Andreev states. Therefore, the measured electrical signature such as the conductance, may not be the manifestation of MBSs. It is important to distinguish the signature of MBSs from the weak coupling between the tunnel probe at the end of the wire and trivial Andreev states in the gated region of the wire.

Zoom In Zoom Out Reset image size

A recent experimental paper by Albrecht et al [30] reports the apparent observation of exponential topological protection in Majorana nanowires. The observation can be explained by the overlap of the two MBSs localized at the two ends of the nanowire is exponentially small as the increase of the wire length. However, the magnetic field dependence of the MBS energy splitting indicates a damping oscillation, which contradicts with the theoretical predictions. To explain this disagreement, some possible scenarios are proposed [42–45]. One of them suggested that for an inhomogeneous nanowire [43–45], the appearance of the trivial Andreev bound states plays an essential role. Our results indicate that in a simple scenario, the trivial Andreev states are hard to survive in the case of a high potential barrier imposed at one end of the nanowire, which is previously expected to prevent the leakage of MBSs. The wave function of MBSs is exponentially suppressed in the barrier region and is hard to leak out of the topological region. However, there always exists a near-zero-energy Andreev bound states in the gated well region. The measured conductance signature about MBSs are unavoidable to coalesce with the feature of Andreev bound states. It is interesting that in a deep gated well region, both the spin-up and spin-down components of the trivial Andreev bound states make considerable contribution in the wave function distribution. In contrast, the MBSs are almost fully spin-polarized. Such a feature could be detected experimentally [26, 46].

In summary, we present an exact solution of the continuum BdG equation for MBSs generated in a superconductor–semiconductor hybird topological nanowire using the trial wave function method. In the simple homogeneous case, it is found that the zero-energy mode could also be produced in the topological trivial regime for a shorter nanowire, which also have oscillatory dependence on the chemical potential. The behaviors of Majorana energy splitting could be well explained by an analytic formula we derived using the trial wave function method.

In the inhomogeneous case, we investigate the low-energy spectra as a function of the chemical potential of the gated region and probability density distribution along the nanowire. Our numerical results show that a chemical potential barrier can prevent the leakage of MBSs and suppress the appearance of extra near-zero-energy Andreev bound states. In contrast, for a chemical potential well case, the trivial Andreev bound states are always generated and coalesced with MBSs at the interface of two parts. Different from MBSs which are almost fully spin-polarized, both spin components of Andreev bound states are significant in the well region. Such a feature could be served as a probe to discriminate the MBSs and Andreev bound states.

We acknowledge the helpful discussion with Prof Shun-Qing Shen and Prof Hai-Zhou Lu. This work was supported by the Natural Science Foundation of China under Grant No. 61474018 and No. 11574173, the Fundamental Research Funds for the Central Universities (No. ZYGX2019J100) and the Open Project of State Key Laboratory of Low-Dimensional Quantum Physics (Grant No. KF202008).

In this appendix, we present a detailed derivation of equation (6). For any z_n, there exists another root −z_n and its wave function is given by

$$\begin{equation}{u}_{n}\left(x\right)={\rho }_{n}\enspace {\mathrm{e}}^{-{z}_{n}x}=\left(\begin{pmatrix}\hfill -{u}_{n1}\hfill \\ \hfill {u}_{n2}\hfill \\ \hfill -{u}_{n3}\hfill \\ \hfill {u}_{n4}\hfill \end{pmatrix}\right){\mathrm{e}}^{-{z}_{n}x}.\end{equation} \tag{ A.1 }$$

Therefore, the total wave function can be expressed as

$$\begin{align}\hfill {\Psi}\left(x\right)& =\sum _{n=1}^{4}\left({a}_{n+}{\rho }_{n+}{\mathrm{e}}^{{z}_{n}x}+{a}_{n-}{\rho }_{n-}{\mathrm{e}}^{-{z}_{n}x}\right)\hfill \\ \hfill & =\sum _{n=1}^{4}\left(\begin{pmatrix}\hfill \left({c}_{n+}{\mathrm{e}}^{{z}_{n}x}-{c}_{n-}{\mathrm{e}}^{-{z}_{n}x}\right){u}_{n1}\hfill \\ \hfill \left({c}_{n+}{\mathrm{e}}^{{z}_{n}x}+{c}_{n-}{\mathrm{e}}^{-{z}_{n}x}\right){u}_{n2}\hfill \\ \hfill \left({c}_{n+}{\mathrm{e}}^{{z}_{n}x}-{c}_{n-}{\mathrm{e}}^{-{z}_{n}x}\right){u}_{n3}\hfill \\ \hfill \left({c}_{n+}{\mathrm{e}}^{{z}_{n}x}+{c}_{n-}{\mathrm{e}}^{-{z}_{n}x}\right){u}_{n4}\hfill \end{pmatrix}\right),\hfill \end{align} \tag{ A.2 }$$

where z_n (n = 1, 2, 3, 4) are the roots with Re[z_n] ⩾ 0. Defining the functions

$$\begin{equation}{S}_{{\pm}}\left(x\right)=\frac{{\mathrm{e}}^{{z}_{n}x}{\pm}{\mathrm{e}}^{-{z}_{n}x}}{{\mathrm{e}}^{{z}_{n}L/2}{\pm}{\mathrm{e}}^{-{z}_{n}L/2}},\end{equation} \tag{ A.3 }$$

one can get

$$\begin{equation}{\mathrm{e}}^{{\pm}{z}_{n}x}=\mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{S}_{+}\left(x\right){\pm}\mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{S}_{-}\left(x\right).\end{equation} \tag{ A.4 }$$

At the boundary x = ±L/2, the wave function equals to zero Ψ(±L/2) = 0. Because S_±(L/2) = 1 and S_±(−L/2) = ±1, we get the boundary condition

$$\begin{equation*}\left(\begin{pmatrix}\hfill \sum _{n=1}^{4}\left({\pm}\mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n1}{d}_{n+}+\mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n1}{d}_{n-}\right)\hfill \\ \hfill \sum _{n=1}^{4}\left(\mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n2}{d}_{n+}{\pm}\mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n2}{d}_{n-}\right)\hfill \\ \hfill \sum _{n=1}^{4}\left({\pm}\mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n3}{d}_{n+}+\mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n3}{d}_{n-}\right)\hfill \\ \hfill \sum _{n=1}^{4}\left(\mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n4}{d}_{n+}{\pm}\mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n4}{d}_{n-}\right)\hfill \end{pmatrix}\right)=0\end{equation*}$$

with d_n± = (a_n+ ± a_n−).

From above boundary condition, we have

$$\begin{equation}\left(\begin{pmatrix}\hfill \sum _{n=1}^{4}\enspace \mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n1}{d}_{n+}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n2}{d}_{n+}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n3}{d}_{n+}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n4}{d}_{n+}\hfill \end{pmatrix}\right)=0\end{equation} \tag{ A.5 }$$

and

$$\begin{equation}\left(\begin{pmatrix}\hfill \sum _{n=1}^{4}\enspace \mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n1}{d}_{n-}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n2}{d}_{n-}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{cosh}\enspace \frac{{z}_{n}L}{2}{u}_{n3}{d}_{n-}\hfill \\ \hfill \sum _{n=1}^{4}\enspace \mathrm{sinh}\enspace \frac{{z}_{n}L}{2}{u}_{n4}{d}_{n-}\hfill \end{pmatrix}\right)=0.\end{equation} \tag{ A.6 }$$

To obtain the nonzero solution of d_n±, the secular equations about E becomes

$$\begin{equation}\text{Det}\left({\mathbf{M}}^{+}\right)\cdot \text{Det}\left({\mathbf{M}}^{-}\right)=0,\end{equation} \tag{ A.7 }$$

where

$$\begin{equation*}{\mathbf{M}}^{{\pm}}=\left(\begin{pmatrix}\hfill {u}_{11}{t}_{1}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{21}{t}_{2}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{31}{t}_{3}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{41}{t}_{4}^{{\pm}}\hfill \\ \hfill {u}_{12}\hfill & \hfill \hfill & \hfill {u}_{22}\hfill & \hfill \hfill & \hfill {u}_{32}\hfill & \hfill \hfill & \hfill {u}_{42}\hfill \\ \hfill {u}_{13}{t}_{1}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{23}{t}_{2}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{33}{t}_{3}^{{\pm}}\hfill & \hfill \hfill & \hfill {u}_{43}{t}_{4}^{{\pm}}\hfill \\ \hfill {u}_{14}\hfill & \hfill \hfill & \hfill {u}_{24}\hfill & \hfill \hfill & \hfill {u}_{34}\hfill & \hfill \hfill & \hfill {u}_{44}\hfill \end{pmatrix}\right)\end{equation*}$$

with ${t}_{m}^{+}=\mathrm{tanh}\left({z}_{m}L/2\right)$ and ${t}_{m}^{-}=c\enspace \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left({z}_{m}L/2\right)$.