Assessing bound states in a one-dimensional topological superconductor: Majorana versus Tamm

Majorana bound states in topological superconductors have attracted intense research activity in view of applications in topological quantum computation. However, they are not the only example of topological bound states that can occur in such systems. We here study a model in which both Majorana and Tamm bound states compete. We show both numerically and analytically that, surprisingly, the Tamm state remains partially localized even when the spectrum becomes gapless. Despite this fact, we demonstrate that the Majorana polarization shows a clear transition between the two regimes.


Introduction
Non abelian quasiparticles have attracted intense theoretical and experimental research activities in the last years, in view of their potential applications in topologically protected quantum computation [1]. Several proposals for the experimental realization or detection of such exotic states of matter in interacting systems have hence been put forward [2][3][4][5][6][7][8][9][10][11][12][13]. Despite steady progress, the need to enforce the coexistence of superconductivity and electronic interactions makes their experimental realization challenging [14][15][16]. In the absence of interactions, the possibility to engineer non abelian Majorana zero modes [17] in spinless p-wave superconductors has been predicted [18,19]. The unambiguous experimental realization of such states is again far from straightforward, due to the fact that spinless p-wave superconductivity is very rare in nature. However, it has been realized that this kind of pairing could be achieved by means of proximity induced superconductivity on semiconducting structures with pronounced spin dependent properties. One of the first example of a system that, once proximitized, can host Majorana bound states, is represented by topological insulators [20][21][22], where strong spin-orbit coupling gives rise to inverted bands behaviour. However, it additionally requires the presence of magnetic barriers to allow for Majorana bound state formation. The difficulty in implanting magnetic barriers has however slowed down the development of such a platform. Recently, the possibility to employ quantum point contacts [23] instead of magnetic barriers has been considered, providing a promising alternative for topological insulator-based Majorana platforms [24,25]. Several other platforms have been proposed, ranging from ferromagnetic atomic chains with superconducting pairing [26], planar Josephson junctions in 2D semiconductors [27][28][29], and the more popular setup based on spin-orbit coupled quantum wire with proximity superconductivity and applied magnetic fields [30,31].
In the last case, several experiments have reported signatures compatible with the presence of Majorana zero modes [32]. However, a clear evidence of the quantization of the zero bias peak, which is one of the predicted signatures ascribed to the presence of Majorana zero modes, or other unambiguous fingerprints of their formation are still lacking to date. Indeed, an additional challenge in the field is posed by the fact that Majorana zero modes behave, in many terms, similarly to more conventional Andreev bound states, which do not exhibit non-abelian statistics, and are thus less useful in quantum computation [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. Indeed, on one hand it is still possible to distinguish between trivial Andreev bound states that emerge at finite chemical potential in clean systems and Majorana bound states due to the fact that the phase containing Majorana fermions and the one with the trivial bound states are separated by a region with extended states only [51]. On the other hand, the distinction becomes more subtle in the presence of disorder, spatial variations of the confinement potentials, or whenever the formation of quantum dots [52] takes place. A detailed analysis of the scenarios in which Majorana fermions and other bound states coexist is hence in order.
In this work, we concentrate on this aspect by studying a simple model where different bound states may coexist. To this end, we consider a one-dimensional finite size spinless p-wave superconductor in the presence of a spatially-periodic modulation of the local chemical potential. The system is numerically investigated by means of extensive exact diagonalization. We find that the bulk of the system is always gapped at the chemical potential, except for a gapless point that marks the boundary between two phases dominated either by the periodic local potential (A) or the superconducting one (B). In the A phase either a bound state is localized at one end of the system or no boundary states are present, depending on the phase of the modulated potential. The possible boundary state is adiabatically connected with the system in the absence of superconductivity and is hence qualitatively speaking a Tamm state [53][54][55]. In the B phase, on the other hand, one recovers the more conventional Majorana scenario, with the zero mode split into two Majorana modes at both ends of the system. Surprisingly, at the transition between A and B the Tamm state remains partially localized. This behavior poses the question if the Majorana fermions in the B phase are still topological in nature or if the reminiscence of the Tamm state eventually spoils their properties. In the second part of the work we answer this question, deriving and analyzing a linearized low-energy model which can be analytically solved. This allows us to interpret and discuss the physics observed in the first part in terms of Goldstone-Wilczek charges [54,[56][57][58][59][60][61][62][63]. Finally, to elucidate the interplay between the different kinds of bound states, and to possibly discriminate between them, we evaluate the Majorana polarization [64][65][66] associated to the boundary states. This quantity is one of the standard tools introduced to characterize topological superconductors. This allows us to conclude that the phase B is indeed characterized by the presence of Majorana fermions.
In more details, the article is divided as follows. In Sec. 2 we analyze the results for the p-wave superconductor model, in Sec. 3 we derive the exactly solvable model and discuss it. Sec. 4 contains the calculation of the Majorana polarization. Our conclusions are finally drawn in Sec. 5.

The quadratic model 2.1. Hamiltonian
The system we inspect is a finite-size, spinless one-dimensional p-wave superconductor in the presence of an additional periodic potential [67]. More specifically the Hamiltonian, defined on a segment of length L, is with the Bogoliubov-de Gennes (BdG) Hamiltonian density H(x) given by (h = 1) and Ψ(x) = (ψ(x), ψ † (x)) T , with ψ(x) a fermionic annihilation field. We impose open boundary conditions, implying Ψ(0) = Ψ(L) = 0. In the Hamiltonian, m is the effective mass, µ 0 sets the filling with 1 k F = 2mµ 0 , ∆ the strength of the p-wave superconducting pairing and τ i (i = 1, 2, 3) are the Pauli matrices. For future use we also introduce τ 0 as the 2X2 identity matrix. Finally we set with the phase ϕ kept as a free parameter. Here, µ 1 (x) is a periodic potential, that can emerge due to external gates [68] or via a coupling to the phonons [69]. The parameter B parametrizes the strength of the term, while the phase ϕ, which is clearly an irrelevant parameter when periodic boundary conditions are imposed, becomes essential in the case of open boundary conditions. The manipulation of the phase ϕ can be envisioned if µ 1 (x) is due to external finger gates [68]. In this case, the spatial variations of the potential can be fully manipulated.
The energy spectrum and the wavefunctions are evaluated via a numerical exact diagonalization procedure, which consists of expanding the eigenfunctions ψ ν (x) of the Hamiltonian on the basis of the states of a free particle in a infinite 1D box of length L and diagonalizing the associated matrix of the Hamiltonian on this basis. This allows us to obtain the energy spectrum E ν and the weights c n (ν) of the eigenfunctions of the problem. It is here worth to mention that since we adopt the Bogoliubov-de Gennes formalism we artificially double the spectrum by introducing a redundant chiral symmetry, and hence only half of the eigenfunctions of the Hamiltonian matrix have physical meaning. Indeed, by only taking the positive energy eigenstates one recovers the correct excitation spectrum. The diagonalization is performed in Mathematica (TM) and up to 450 box states per BdG sector are employed depending on the values of B and ∆, to ensure numerical convergence with a relative error δ 10 −3 on the energy spectrum.

Tamm states
We begin our analysis by considering the case ∆ = 0. Figure 1 shows the energy spectrum for different values of the strength B and different phases ϕ. As a general feature we observe that for B = 0 the periodic potential opens a gap at zero energy -where the chemical potential is set. Also, secondary gaps can occur when the amplitude B becomes large, B µ 0 . This additional gap, visible in 1(a) for ν ∼ 150, opens, having in mind periodic boundary conditions and a folded scheme for taking into account the periodic perturbation, at k = 0. For the same magnitudes of B, an additional feature appears for small ν (see the blue and red curves). This is due to the change in nature of the lowest energy states, which, for B > µ 0 , are dominated by the periodic potential. Comparing Figs.1(a) and (b) one can observe that for almost all the states, the overall shape of the energy spectrum depends weakly on the phase ϕ, being almost insensitive to it for E µ 0 . However, inspecting the case ϕ = π/2 one can observe that a zero-energy mode occurs within the gap, which suggests that the low-energy sector of the spectrum may exhibit a more pronounced dependence on ϕ. To check this fact, Fig. 2(a) shows the behaviour of the first eight eigenstates with E ν ≥ 0 as a function of ϕ for the representative case of B = µ 0 /2. As one can clearly see, the lowest energy state exhibits a quite dramatic dependence on the phase, becoming a zero mode for ϕ = (2p + 1)π/2 with p an integer. The wavefunction of such a zero-energy state is shown, for ϕ = π/2 in Fig. 2(b) and it corresponds to a localized state, which in this context is usually referred to as a Tamm state. Such state is respectively localized either at the right (ϕ = π/2 + 2π p) or at the left (ϕ = 3π/2 + 2π p) edge of the system. Such Tamm state is starkly different from the delocalized states obtained at ϕ = π p, which extend over the entire length of the system -see the representative example reported in Fig. 2(c) for the case ϕ = 0. In general, for ϕ = π p a bound state within the gap always occur, with possibly a non-zero energy and a larger localization length.

Majorana states
The regime B = 0 with ∆ = 0 is more known since it represents the paradigmatic model that hosts Majorana bound states: even in this case a gap opens around the chemical potential and a fermionic zero-energy mode appears in the spectrum. The BdG wave function corresponding to this fermionic state is significantly non-zero close to both ends of the segment -see Fig. 3. Formally, one can interpret this non-local fermionic state as two local Majorana states, each one localized close to one end only. This is obviously meaningful only when the localization length is smaller than the length of the system.

Competition between Tamm and Majorana states
Let us now address the case when both B and ∆ are non-zero and thus the two gaps induced by the oscillating potential and by the p-wave pairing compete. We are particularly interested into the interplay between the two (qualitatively different) localized Tamm and Majorana states and thus we set from now on ϕ = π/2 to achieve the most dramatic effects. Indeed, ϕ = π/2 corresponds to a maximally localized Tamm state. As a general feature, we find that in this case, for any value of B, ∆ a zero-energy mode occurs.    Figure 4(a-e) shows the probability amplitude |ψ 0 (x)| 2 of this zero-energy mode for B = µ 0 /4 and increasing values of the superconducting pairing. As ∆ is increased from zero, the localized state broadens -see panel (b). However, when ∆ = B -see panel (c) -a peculiar phenomenon occurs: the former Tamm state gets once again sharper while the wavefunction in the rest of the system becomes non-zero and with a flat envelope. When ∆ > B, see panel (d), a second localized state develops. Finally, when ∆ B the wavefunction tends to recover a symmetric shape localized at both ends of the system. To better understand what is going on, we study the spectral weight w(n, E) = |c n [ν(E)]| 2 , shown as a density plot as a function of the box state index n and the energy E in Fig. 5 2 for energies around E = 0. In general, this tool allows us to get a glimpse of the structure of the energy spectrum of the system. Indeed, as a warm up panel (a) shows the case B = ∆ = 0: as one can see, the spectral weight is sharply peaked on one box state per energy value and the yellow traces faithfully reproduce the energy spectra of the particle (E > 0) and hole (E < 0) sectors of the BdG picture for the Hamiltonian shows the situation for B = µ 0 /4 and ∆ = 0, when a Tamm state occurs: The energy gap is clearly visible and the blurry feature at E = 0 is the spectral representation of the Tamm state, whose amplitude is smeared across many delocalized box states in order to produce a state localized in space. The situation is somewhat similar for ∆ > B -see panel (d) -again consistent with the picture of a wavefunction with sharp peaks at the system ends. Very peculiar is however the situation for B = ∆, shown in panel (c): at this transition point the spectrum is gapless and exhibits two linear modes around E = 0. Thus, the zero-energy mode is embedded into a gapless spectrum and yet it still represents a bound state. To summarize our findings so far, even though a gap closing is present at the level of the bulk spectrum, bound states can exist for every value of B and ∆. However, assessing the nature of such bound states (Tamm-or Majorana-like) is not a trivial affair in the context of the full quadratic model which can only be solved numerically. A notable exception is represented by the point ϕ = 0. In that case, indeed, no Tamm state is present and a bound state is necessarily a Majorana state. However, the general case needs further inspection. In order to tackle this problem and provide an answer, we now proceed developing an analytically solvable low-energy model which is however able to capture all the features described so far.

Linearization
To develop an effective, low-energy model valid around the chemical potential we consider a linear dispersion relation instead of the quadratic one. This is meaningful as long as µ 0 is the largest energy scale involved [70][71][72]. In details, we approximate the fermionic operator as [73][74][75][76] ψ with ψ R,L (x) fermionic operators for right-or left-moving electrons. In particular, by using the particle in a box states one can write and ψ L (x) = −ψ R (−x) ensuring, together with the 2L periodicity of the fields, that open boundary conditions ψ(0) = ψ(L) = 0 are satisfied. Here, the fermionic operators d n are associated to the n−th eigenstate of a particle in a box. Note that, as an approximation, the sum is here extended to all integers and, consequently, a slight notation abuse has been adopted having defined the operators d n for all the integer values n. This approximation scheme is always performed when discussing the bosonization procedure, and is discussed, for example, in Ref. [81]. The kinetic energy is then approximated with K L with v F = (πn F )/(mL) representing the Fermi velocity. In order to write the full Hamiltonian within the linear approximation, it is useful to introduce the enlarged BdG spinor By imposing that µ 0 is the largest energy scale (we indeed neglect terms proportional to 1/L with respect to terms proportional to k F ), we finally get the approximate form for the Hamiltonian with In order to interpret the boundary conditions and the peculiar effects they may have on the bound states, it is useful to think in the following way. Consider a pair of chiral fermions χ R (x) and χ L (x), independent from each other and organized in the BdG spinor

defined on the whole real axis, and with Hamiltonian
where we have introduced a localized backscattering potential with m 0 the backscattering strength. It can be easily shown that in the limit m 0 /v F → ∞ the wavefunctions in the region 0 < x < L become disconnected and acquire the boundary conditions χ L (0) = −χ R (0), χ L (L) = −χ R (−L) [77][78][79]. The Schroedinger problem associated to H χ , when considered for 0 < x < L and when m 0 /v F → ∞ is hence perfectly equivalent to the one related to H in Eq. 9. The advantage in resorting to H χ is conceptual, since it allows to give a simple interpretation to the results described in the previous section. Indeed, from a physical point of view, H χ is equivalent to the Hamiltonian of a quantum spin Hall liquid gapped by both superconductivity, with gap ∆, and a magnetic mass pointing in the direction given by φ in the presence of two strong backscattering centers with opposite magnetization localized at x = 0, L.

Qualitative interpretation of the results
The first insight into the results comes from the computation of the spectrum as obtained by neglecting the presence of boundaries. If we consider the Hamiltonian H ∞ given by where the operators are now being defined, with a slight abuse of notation, on the whole real axis, one can obtain the energy spectrum. The spectrum is promptly found to be , 2} and is shown in Fig. 6. As can be seen, the spectrum is always gapped for all B > 0, ∆ > 0 except for B = ∆ -see panel (b)where 1 (k) is gapless. The phase space point B = ∆ marks a quantum phase transition where the gap changes from magnetic (B > ∆) to superconducting (B < ∆). It is however crucial to notice the presence of a secondary gap of magnitude B + ∆ in 2 (k) which is always present for every B > 0, ∆ > 0. Notice also the qualitative similarity between these spectra and the spectral features found in the quadratic model substantiating the validity of the linear approximation at low energies, compared to µ 0 .
Turning to the inspection of bound states it is crucial to observe that the addition of a confining term proportional to m 0 produces a magnetic gap in the spectrum. It is well known [80] that in the superconducting dominated regime (B < ∆) Majorana bound states appear at the boundaries of the system. However, in the opposite case, when B > ∆ the situation is more subtle, as this model would implement the physics of so called Goldstone-Wilczek fractional solitons [62]. Since this regime is adiabatically connected to ∆ = 0, in order to build some physical intuition one can refer to the case in which superconductivity is absent. The system is then equivalent to a quantum spin Hall liquid, in the absence of superconductivity, gapped by two magnets of equal strength at x = 0, L and characterized by different magnetization angles: it hence hosts bound states. Assuming the decay length of these bound states to be shorter than the system size L, the two bound states are essentially independent, and thus it is sufficient to concentrate on a single boundary to capture the relevant physics.

Linear model with one boundary and ∆ = 0
We thus inspect a simpler linear model in the regime ∆ = 0, L → ∞ considering the bound state localized around x = 0. Its wavefunction ψ m (x), that since superconductivity is now absent (∆ = 0) has only the two right-and left-moving components, can be easily obtained and reads The solution is only acceptable (normalizable) for π < φ < 2π. The corresponding energy is 0 = −m cos(φ). Note that for φ = 3π/2 the state, mimiking the Tamm state of the quadratic model, is at zero energy. The scenario just derived perfectly agrees with the behavior of the left boundary of the quadratic model, although in that case the wavefunction in Eq. 14 by means of Eq. 5 has one component only, and is rebuilt from the condition in Eq. 5. Moreover, in the quadratic model, the unphysical particle-hole symmetry of the Bogoliubov-de Gennes equation is present while here the redundancy [82,83] is not needed since superconductivity is not present. We hence have a single bound state instead of two at opposite energies. The addition of a second physical boundary at x = L produces effects that can be now easily understood by leveraging on the single boundary case just discussed. Indeed, the bound state at the boundary at x = L is present when sin(φ) > 0 -that is for 0 < φ < π -when the bound state is absent at x = 0. Such exponentially increasing solution is admissible due to the added barrier at x = L -which prevents the divergence of the wavefunction -and its decay through the dot. This again is perfectly compatible with the results of the quadratic model discussed in Sec. 2. A final comment for this section is the behavior at φ = 0, π. For those values of the angle, the energy of the bound state reaches the bulk bands and the bound state becomes delocalized, and hence, effectively, disappears. In the following we will concentrate on the case in which the bound state is maximally localized around x = 0, that is for φ = 3π/2.
3.4. The wavefunction of the bound state for ϕ = 3π/2 A further interesting behavior of the bound state found within the quadratic model is the fact that, when present, it does not become completely delocalized at the quantum phase transition B = ∆. To confirm that this fact occurs also within the linearized model, we here report the BdG expression of the wavefunction ψ 0 (x) of the bound state characterizing the linearized model in Eq. 9, together with the condition in Eq. 5, for φ = 3π/2 and arbitrary B > 0, ∆ > 0. We find where Both the particle and the hole components u(x) and v(x) are characterized by a profile that is peaked around the edge even at the transition point B = ∆. This case is shown in Fig.  7, where a plot of |u(x)| 2 and |v(x)| 2 is reported for B = ∆: one can clearly notice that at the transition point, although the wavefunction is sharply localized at x = 0, it does not completely decay exponentially to zero. Also this fact is in excellent agreement with the results obtained in the full quadratic model. As apparent from the analytical form of the wavefunction, this fact is due to the presence of the B + ∆ energy scale. A natural question hence arises: are the Majorana zero modes, in the superconducting dominated phase B < ∆, influenced by the persistence of the Tamm state across the quantum phase transition?

Majorana polarization
The Majorana polarization P M [64] is a quantity that has been proposed in order to discriminate between Majorana fermions and other types of Andreev bound states. We now evaluate the Majorana polarization associated to the state ψ 0 (x), defined as This absolute value of P M is close to one (in the large gap limit) if the bound state is a Majorana bound state, while it reaches different values for regular Andreev bound states, which allows to assess the character of ψ 0 (x). In the case under examination we find Looking at the previous expressions, two limiting cases can be extracted: for B = 0 one finds P M = tanh L∆ v F , with P M → 1 when ∆ → ∞.  The Majorana polarization is shown in Fig. 8.

Conclusions
In this work we have characterized a simple model for a finite size, one dimensional topological superconductor in the presence of a competing normal gapping mechanism arising from a position dependent potential. The first part of the article deals with the numerical inspection of the model. We find that the model has two phases, one characterized by Majorana fermions at both the ends, and the other by at most a single Tamm bound state localized at one side. The energy of the bound state, when present, strongly depends on the phase of the position dependent potential. Surprisingly, we find that the wavefunction of the bound state remains strongly peaked close to the edges even at the gap closing point separating the two phases. In the second part of the article, we develop an exactly solvable model for interpreting the results of the first part. Finally, in the last section we evaluate the Majorana polarization for a fully developed bound state of the model and show that, despite the fact that the Tamm state is not completely delocalized at the transition point, deep in the superconducting phase the Majorana bound states are not strongly affected by the position-dependent potential.