Zero modes of the Kitaev chain with phase-gradients and longer range couplings

We present an analytical solution for the full spectrum of Kitaev's one-dimensional p-wave superconductor with arbitrary hopping, pairing amplitude and chemical potential in the case of an open chain. We also discuss the structure of the zero-modes in the presence of both phase gradients and next nearest neighbor hopping and pairing terms. As observed by Sticlet et al., one feature of such models is that in a part of the phase diagram, zero-modes are present at one end of the system, while there are none on the other side. We explain the presence of this feature analytically, and show that it requires some fine-tuning of the parameters in the model. Thus as expected, these `one-sided' zero-modes are neither protected by topology, nor by symmetry.


I. INTRODUCTION
One of the characteristic features of many topological phases is the presence of gapless boundary modes. The (fractional) quantum Hall states are a prime example 1-3 , and their boundary modes provide strong evidence of the topological nature of these states. Another prime example is the Kitaev chain, whose topological p-wave superconducting phase features so-called 'Majorana zero modes' at its edges. 4 Trying to establish the existence of the topological phase is often done by trying to establish the presence of the boundary modes. This has led to strong evidence for the topological phase in for instance strongly spin-orbit coupled nano-wires that are proximity coupled to an s-wave superconductor in the presence of a magnetic field [5][6][7][8][9] , or in chains of magnetic ad-atoms [10][11][12][13] . It has been proposed that the zero energy Majorana bound states can be used as topologically protected qbits, for quantum information processing purposes 14,15 . By now, there exist various proposals to manipulate these q-bits, either in T-junction systems, in which the Majorana bound states can be braided explicitly 16 , or in Josephson coupled Kitaev chains, in which the coupling of the various chains allows operation on the q-bits 17 .
Despite the intense research on the Kitaev chain models, there are still interesting features that deserve attention. In this paper, we look into one of them. It was observed by Sticlet et al. 18 , that the zero-modes of Kitaev chains carrying a current, i.e., in the presence of a gradient in the phase of the order parameter, have interesting properties. The most striking feature is that is it possible that at one edge of the chain, there is pair of Majorana bound states (or better, one 'ordinary' Dirac zero mode), while there is no zero mode at the other end of the chain. Clearly, from a topological point of view, this means that the chain is in a trivial phase, but it is nevertheless worthwhile to investigate these zero-modes further. In this paper, we explain the presence of these zero-modes, via an exact solution of the zero modes of an extended Kitaev chain, i.e., in the presence of both complex and next nearest-neighbor hopping an pairing terms. We show that it is necessary to fine tune the couplings in order that these 'one-sided Dirac modes' to exist, but under these fine-tuned conditions, they can only disappear if the bulk gap closes signaling a phase transition or a crossover. Dropping the fine-tuning will gap out these zero modes immediately, leaving behind low-energy subgap modes. Apart from the analytical solution of the zero modes, we also present the solution of the full spectrum of the open Kitaev chain, for real, but otherwise arbitrary parameters, which does not seem to have appeared in the literature before.
The outline of the paper is as follows. We start in Sec. II by a brief review of the Lieb-Schultz-Mattis method to solve open quadratic fermionic systems, and focus on the case of complex couplings, which is essential for our purposes. In Sec. III, we provide the full solution of the open Kitaev chain, with real, but otherwise arbitrary couplings. In Sec. IV, we study the effect of next nearest-neighbor and complex couplings. Here, we focus entirely on the exact solutions for the zero-modes, and start by considering the effects of next nearest-neighbor couplings and complex pairings separately, before coming to the most interesting case, when both are present. In Sec. V, we discuss the results of the paper. Some details are delegated to the appendices.

II. THE LIEB-SCHULTZ-MATTIS METHOD
In this paper, we study the zero modes of Kitaev-like chains in the presence of longer range hopping and pairing terms, specifically next nearest neighbor (NNN) ones. In particular, we are interested in the case where these couplings are complex. To study these systems, we use the method has been introduced by Lieb, Schultz and Mattis (LSM) 19 who used it to solve the XY chain, for various types of boundary conditions. For a quadratic fermionic Hamiltonian with periodic boundary conditions (PBC), one diagonalizes the Hamiltonian by using a Fourier transformation, followed by a Bogoliubov transformation in the case of superconducting model. Without translational invariance one can still perform a Bogoliubov like transformation directly in real space. It was this method that LSM used to find the spectrum of the open XY chain (after using a Jordan-Wigner transformation to transform the spin degrees of freedom to polarized arXiv:1709.00959v2 [cond-mat.mes-hall] 11 Apr 2018 fermions).
In this section we review the LSM method and follow their notation for convenience. We consider two different cases. First, we look at the Hamiltonian with real couplings and recall how one can derive the spectrum of the model analytically. Second, for a general quadratic Hamiltonian with complex couplings we present the equations governing the zero mode solutions, which we use throughout the paper.
Following LSM 19 , we consider the general quadratic Hamiltonian of polarized fermions as follows, in which c i is a fermion annihilation operator on site i, A is a hermitian matrix, B is an antisymmetric matrix and N is the number of sites. Using a Bogoliubov like transformation, one can define new fermion operators, and diagonalize the Hamiltonian: in which α labels the states and g α,i and h α,i are two functions, which are to be determined. This transformation is canonical, in the sense that new operators obey the fermionic anti-commutation relations, i.e. {η α , η † β } = δ αβ .
Using the equation of motion, [H, η α ] = −Λ α η α , one finds the equations for g α,i and h α,i : In order to find the full spectrum of the Hamiltonian, we now consider the case for which A and B have real elements. In this case, one defines new variables as which we combine into row vectors, φ α = (φ α,1 , . . . , φ α,N ) and ψ α = (ψ α,1 , . . . , ψ α,N ). Summing and subtracting Eqs. (4) and (5) gives two coupled equations for φ α and ψ α , We note that the matrices act from the right side on the vectors. By acting with A + B on Eq.(8) and A − B on Eq.(9) from the right, the equations decouple To find all the eigenvalues Λ α and states η α , one solves these two decoupled equations for φ α and ψ α . We explain how to do this in more detail in the next section for the open Kitaev chain 4 with real, but otherwise generic parameters.
It is well-known that fermionic systems can host Majorana zero modes on the edges of the system, which signals that the system is in a topological phase. In this paper, we study the zero modes of Hamiltonians with complex parameters, so we now allow the matrices A and B to be complex again. To distinguish a Majorana mode from the ordinary modes, we use stared labels, such as α * . The Majorana modes satisfy η α * = η † α * . For a finite system, the energy of a Majorana mode is exponentially small in the system size; for instance in the case where we have a system with N sites the energy scales as Λ α * ∼ e −κN with κ > 0 4,19 . Hence we are interested in finding general equations which allows one to find the corresponding states with zero energy, i.e. Λ α * = 0, in the thermodynamic limit.
We thus search for a Majorana solution of Eqs.(4) and (5) with zero energy. Setting h α * ,i = g * α * ,i in Eqs. (4) and (5) gives: By summing and subtracting these equations we get, We use these equations to explore the wave functions (g α * ,i ) of the zero modes in different cases in the following sections. Before closing this section, we write the η operators in terms of Majorana operators for future reference. Using φ and ψ as defined above and defining Majorana operators as γ A,j = c † j + c j and γ B,j = i(c † j − c j ), we write the fermion annihilation operator as follows The algebra of Majorana operators can be calculated from the canonical anti-commutation relations of the c operators, Specifically, for the zero mode solution we can write the corresponding fermionic operator as follows:

III. THE SPECTRUM OF THE OPEN KITAEV CHAIN
In this section, we use the method of LSM to find the full spectrum of the Kitaev chain 4 , for an open chain, with real parameters, in particular we consider Here, µ denotes the chemical potential, ∆ the strength of the pairing term, and we chose the hopping parameter t = −1 1 .
Despite the fact that this model has been studied thoroughly, these results do not seem to have appeared in the literature, and we will use it to set the notation. Because we are interested in the zero-modes of more generic situations in the remainder of the paper, we also quickly review the nature of the zero-modes. These latter results are not new, but appeared in [19][20][21] and for generic parameters recently in 22,23 .
It is well known 4 that the Kitaev chain is in a topological phase for |µ| < |t| and ∆ = 0. A profound feature of topological phase is the presence of a Majorana zero modes, that are exponentially localized near the edges of the system. In addition, the energy associated with this zero mode is exponentially small in the system size.
To set the scene, we follow Kitaev to show the presence of Majorana zero modes, by considering the special case of ∆ = 1 and µ = 0. In this case, the Hamiltonian in terms of Majorana operators becomes, In this Hamiltonian, γ A,1 and γ B,N are absent and therefore commute with it. So one can form a non-local fermionic state, f 0 = 1 2 (γ A,1 + iγ B,N ). The presence of this non-local fermionic mode is the characteristic feature of the topological phase of the Kitaev chain. For ∆ = −1, the unpaired Majorana operators would be γ B,1 and γ A,N , owing to the p-wave nature of pairing.
We leave this fine tuned point and consider arbitrary ∆, but keep µ = 0 for the moment. This corresponds to the XY model, which was solved exactly by LSM for |∆| < 1, that is, the full spectrum including the wave functions were found 19 . For an open chain, there is a state with an exponentially small energy as a function of the system size. The wavefunction of this state is exponentially localized on the edges, namely φ n ∼ 1−|∆| 1+|∆| n where n denotes the position of the site measured from 1 The sign of t is irrelevant for the spectrum, but we set t = −1, because of the simpler relation with the associated XY model as studied in 19 .
the left side of the chain. The associated ψ n is localized on the right edge. Another fine tuned point that was studied previously corresponds to the transverse field Ising model(TFIM), that is t = −1, ∆ = ±1 but arbitrary µ. Pfeuty showed that this model has a Majorana zero mode if |µ| < 1. The associated wave function takes the form φ n ∼ |µ| n and is localized on the left edge of the system 20,21 .
To find the Majorana zero modes for the general case, it is advantageous to first consider the model with periodic boundary conditions. That is, we need to consider the hopping and pairing terms for the last site as well. We denoted the periodic Hamiltonian by H P BC = H + H N where: The solution of the periodic model is well known, and obtained by using a plane-wave ansatz for the wave functions (i.e., by Fourier-transforming the model). Using the method outlined in the previous section, we start by solving Eqs. (10) and (11) to find the spectrum. Since φ and ψ are related via Eqs. (8) and (9), we focus on φ. Writing Eq. (10) gives us one recursion relation: where n denotes the sites and runs from 1 to N . Upon setting φ k,n ∼ e ikn , where we use the momentum k as a label, one finds the eigenvalues: where m runs over 0 to N − 1. If one considers antiperiodic boundary conditions, the dispersion is the same though the allowed values of k change to k = 2π N (m + 1 2 ). To study open chains, we make use the functional form of the dispersion. In addition, by using the LSM method for open chains, we naturally have to consider both the sectors with even and odd number of particles.
We now consider the full spectrum of the open chain. Here, we merely give the results, and refer to App. A, where the details of calculation are presented.
For the open chain we find the same recursion relation in the bulk which is valid for 3 ≤ n ≤ N − 2. However, we also have four boundary equations which should be treated separately (see App. A). We start by dealing with the bulk equations, using the method of LSM. That is, we use the same 'function' for the eigenvalues, though with a generic parameter α instead of the momentum k. To find the allowed values for the parameter α, one uses the 'boundary equations'. Hence we parametrize the eigenvalues as: and α is the label for the state. For the states, we use a power law ansatz, φ α,n ∼ x n α , and we find four solutions, x α = e ±iα and x α = e ±iβ where cos α + cos β = 2µ Note that α and β are not necessarily real, but the way we parametrize x α turns out to be convenient. As described in the App. A, the relevant linear combination that one uses to find a solution for the boundary equations is: in which A 1 and A 2 are constants that are related via Eq. (A24). The boundary equations give rise to another constraint on α and β. This constraint can be shown to take the following form To obtain the full solution of the model, one needs to solve Eqs. (25) and (27) simultaneously. Though this can not be done analytically, it is straightforward to obtain the solutions numerically. Thus, we have characterized all the eigenvalues and eigenvectors φ α,n and by using Eq. (10), one finds ψ α,n . Now we want to study these solutions and see when this model has a Majorana solution and what the corresponding wavefunction is. To find such solutions, we consider thermodynamic limit, i.e. N → ∞, which makes the calculations easier.
We first mention that ∆ can always set to be positive. One way to see this is by considering the transformation under which c j maps to e i π 2 c j . This transformation changes neither the hopping nor chemical potential term, but ∆ changes to −∆. In addition solutions for µ < 0 can be constructed from the solutions for µ > 0. One can take the solution for µ > 0, say (α, β) = (r, s). Now consider (α, β) = (r + π, s + π). This change gives a minus sign for the LHS of Eq. (25) as required, however it leaves Eq. (27) unchanged. Therefore, we restrict ourselves to ∆, µ > 0.
First we look at the solutions for large values of µ. In this case one can see that Eqs. (25) and (27) have N distinct real solutions for α, where we restrict α to lie in the range 0 < α ≤ π (α = 0 gives φ n = 0; for more details, see App. A). However by decreasing chemical potential solution with the smallest value of α will 'disappear'. It is well known that for µ < 1 one real solution is lost in the thermodynamic limit. For a finite chain this happens is a finite size correction. Thus, for µ < 1 one must find an additional, complex solution. To find this solution, we consider three different cases.
1) ∆ < 1 and √ 1 − ∆ 2 < µ < 1: In the thermodynamic limit one can check that the following solution satisfies Eqs. (25) and (27), Furthermore it is straightforward to check that Eq. (24) gives Λ α * = 0, hence the solution is indeed a zero mode. The wave function φ α * ,n for this zero mode is where C is a normalization constant. Moreover, it can be shown that based on structure of A − B and A + B matrices, one has ψ α * ,n = φ α * ,N +1−n . From the fact that ξ 1 < ξ 2 , it follows that φ α * is localized on the left edge while ψ α * is localized on right edge of the system. Hence we found two Majorana operators, that are located on the edges of the system, and the associated wavefunctions decay exponentially.
2) ∆ < 1 and µ < √ 1 − ∆ 2 : In this range, one needs to use a different parametrization if one wants to use real parameters, as is evident from Eq. (29). This parametrization reads Basically we changed one of the characteristic length scales to become a wave vector. As in the previous case, this solution is indeed a zero mode, i.e. Λ α * = 0, whose wavefunction is given by: This result indicates that φ (ψ) is localized on the left (right) edge with an oscillatory decaying wave function. We should point out that this result was obtained earlier by 22 . In addition, it was observed that the correlation functions in the model with PBC are oscillatory in the same regime, i.e., for µ < √ 1 − ∆ 2 with ∆ < 1, see for instance Refs. 24-26.
3) ∆ > 1: For this regime √ 1 − ∆ 2 is imaginary, hence the previous solutions are not applicable. The new root can be written as One can check that this solution represents a zero mode with the wave function: For this solution ξ 1 < ξ 2 since µ < 1 and this guarantees that φ (ψ) is localized on the left (right) edge.

IV. ZERO-MODES FOR NEXT NEAREST-NEIGHBOUR AND COMPLEX COUPLINGS
In this section we study the zero modes in the presence complex hopping and pairing terms, both in the case with nearest neighbor hopping and pairing terms, as well as next-nearest neighbor (NNN) hopping and paring terms. The complex amplitudes model the presence of a phase gradient in the system.
In their fermionic incarnation, these generalized Kitaev models were studied in Refs. 18,27-29. In the language of spin models, adding NNN terms gives rise to so-called (one-dimensional) cluster models 27,30-33 , but we concentrate on the fermionic version of these models.
An important feature of these models is the possibility of having more than one zero modes at each end, which is possible due to the presence of longer range terms. This can also be understood in terms of the classification of topological insulators and superconductors 34,35 . The Kitaev chain with real coupling constants belongs to class BDI for which the different topological phases can be labeled by the elements of Z, in the absence of interactions. Adding interaction changes this picture such that new classification is given by Z 8 instead 36 . In the case with only nearest neighbor hopping and pairing terms, the model describes phases with at most one Majorana mode at each end of the system. However by adding NNN terms one finds phases with two Majorana modes at each end. This means that there would be two distinct topological phases with one and two zero modes solutions (in addition to the trivial phase, which does not have a zero mode).
Proposals for using the non-local fermionic state as a qubit, requires the ability to move Majorana edge states and even to do braiding. One proposal to achieve this is by inducing a phase gradient in the superconductor order parameter, i.e ∆ j = ∆e iθj , with non-uniform θ j 37 . Having a complex superconductor order parameter breaks the time reversal symmetry in which case the model belongs to class D. For class D, we have the Z 2 classification which means that the system could be either in the topological phase with at most one Majorana zero mode at each end, or in the trivial phase. Surprisingly, Sticlet et al. showed that NNN terms with a phase gradient can exhibit an exponentially localized fermionic zero mode on just one edge 18 . Such a phase, though it is not topologically protected, has local zero modes. In Ref. 18 these models were investigated numerically. Here we present an analytical solution and study the zero-modes in detail. We first review the Kitaev chain with NNN terms. After that, we study the effect of a constant phase gradient in the Kitaev chain. Finally, we combine the two complications and consider NNN terms and a phase gradient simultaneously.

A. Next nearest-neighbor couplings
In this section we consider the Kitaev chain and add NNN hopping and pairing terms. We start with the case for which all the parameters are real, hence the Hamiltonian belongs to class BDI. Setting µ = 1, the problem has four energy scales, corresponding to the two hopping and two pairing amplitudes. To simplify the calculation we set the NN hopping and pairing terms equal to each other and we do the same for the NNN terms. Sticlet et al. studied this model under the same assumptions 18 . We consider the model with arbitrary complex parameters in Sec. IV D.
Thus, the Hamiltonian reads, where λ is the NNN hopping and pairing amplitude. To obtain the phase diagram, we first consider the model with periodic boundary conditions 18,27 . We do a Fourier where the τ α are Pauli matrices that act in the Nambu space Ψ k . The Hamiltonian can be written as H k = h(k) · τ . One can find the phase diagram by calculating the winding number for h(k) 18,35 or by looking at gap closing lines 27 . The phase diagram is presented in Fig. 1.
The gap closes along the lines λ = µ + t, λ = µ − t and λ = −µ for |t| < 2|µ|. Note that in the figure we used µ = 1. Before looking at the zero mode solution(s) of an open chain, we first consider some limiting cases to understand the phase diagram. For very small |t|, |λ| |µ| we get the trivial phase. The "0" in Fig. 1 indicates that there are no Majorana zero modes in this part of the phase diagram. Outside of the trivial region on the vertical axis where t = 0 we have two decoupled Kitaev chains, hence there are two zero modes at each end. For a fixed λ, adding NN terms couples these two chains. The two zero modes survive until the gap closes, thereafter there will only be one zero-mode at each end. The horizontal axis with λ = 0 (i.e., the original Kitaev chain) belongs to this later region which is indicated by "1" in the Fig. 1. To find the wave functions of the zero modes, we use Eqs. (8) and (9) with Λ α = 0. From Eq. (16) we see that if η α is a Majorana mode (i.e., η † α = η α ), ψ has to be purely imaginary. So for convenience we define ψ = iψ and we get g = 1 2 (φ + iψ). We use this convention from now on. We obtain the following 'bulk' equations −µφ n + tφ n+1 + λφ n+2 = 0 (39) The 'boundary' equations are , which gives the result where L ± and R ± are real normalization constants (the subscript "0" in the length scales indicates that we deal with a zero phase gradient). We can extract the phase diagram from this result 27 and we set µ = 1 to be able to compare with Fig. 1. For regions where λ > 1 + |t| or both λ < 1 − |t| and λ < −1 (corresponding to the upper and lower regions of Fig. 1), one can see that |x 0,± | < 1. This means that in these regions that are indicated by "2" the system has two independent zero mode solutions. In the right part of the phase diagram where 1 − t < λ < 1 + t, there exists only one zero mode since |x 0,+ | < 1 and |x 0,− | > 1. If 1 + t < λ < 1 − t we also have one zero mode, however, in this case |x 0,− | < 1 and |x 0,+ | > 1. We note that in these regions, the boundary equations are also satisfied in the large N limit.
It is also interesting to note that for t 2 + 4λ > 0 the roots are real. Still they could be negative in some regions which gives rise to an oscillatory behavior of the wave functions, which are then proportional to (−1) n . For t 2 + 4λ < 0 the roots become complex. The red, dashed lines in Fig. 1 specify the upper boundaries of this region (in the case λ < −1). In this case |x ± | = 1 √ |λ| which gives us the the criterion λ < −1 in order to have a zero mode (in the region t 2 + 4λ < 0). In this part of the phase diagram the correlation length only depends on λ, while the NN coupling t only affects the oscillatory part of the wave function. Before moving to the case with both NNN terms as well as with a phase gradient, we first study the Kitaev chain with just a constant phase gradient.

B. Phase gradient in the order parameter
In this subsection, we consider the Kitaev chain, but with a phase gradient in the superconducting order parameter. In the case of a superconductor with a super current, the pairing term has a site dependent phase ∆ j = e ij∇θ where ∇θ is the constant phase gradient, while j indicates the position of the site. In this case, the Hamiltonian reads (42) This Hamiltonian belongs to class D. As we indicated above, the topological phases are labeled by elements of Z 2 , which means that the system could be in the topological phase with one Majorana zero mode at each end. Changing the gauge, we transform the fermionic operators as c j → e ij ∇θ 2 c j . This transformation gives us siteindependent couplings, but now also the hopping parameter has become complex. The transformed Hamiltonian is To find a zero mode solution we use Eqs. (14) and (15).
where L is the normalization constant to make γ 2 L = 1 and the sum is over the sites. This Majorana mode is located at the left side of the system, and is a solution in the large N limit. The right Majorana mode is more complicated, where R is the normalization constant to make γ 2 R = 1 and the sum is over the sites. Using the Majorana modes γ L and γ R , one can construct one fermionic mode f 0 = 1/2(γ L + iγ R ) as usual. We note that to have a localized zero mode we have the criteria |µ| < | cos( ∇θ 2 )|. This means that turing on the phase gradient shrinks the topological region. Second, we see that the left Majorana consists only of γ A Majorana operators (recall the definition above Eq. (16)), however, the right one involves both γ B 's as well as γ A 's. In the case that ∇θ = 0 the left Majorana mode only involves γ A operators and the right Majorana modes only γ B operators. This feature of the solution comes from the fact that for real A and B matrices (see Eqs. (14) and (15)), the equations governing φ andψ are decoupled -recall that g = 1 2 (φ + iψ). Adding the phase gradient makes these matrices complex, hence the equations become coupled and the solutions become more complicated. The direction of the phase gradient shows itself in the elements of the A and B matrices and gives rise to this asymmetry; the "left-right" symmetry is broken explicitly.
In the next section we add NNN terms to the current problem 18 . The results presented in the current subsection are useful to understand zero mode solution(s) when one adds the NNN terms.
C. Next nearest neighbor terms along with a phase gradient in the order parameter We now consider NNN terms in the presence of a constant phase gradient. Again we set the hopping and pairing amplitudes equal to each other for both the nearest neighbors and NNN terms. Following Sticlet et al 18 , the Hamiltonian reads, where we assumed the same phase dependence for the nearest neighbor and NNN pairing terms, with the same phase for both terms involving the first site. As we mentioned above, for ∇θ = 0 this model has a trivial phase without any zero modes and two topological phases that hosts one or two Majorana zero modes respectively (see Fig. 1). For ∇θ = 0, the model belongs to class D. This means that, contrary to ∇θ = 0 case, there is only one type of topological phase. The phase that had two Majorana zero modes becomes trivial upon adding the phase gradient. The natural question is then what happens to the phases with two Majorana edge states? Despite the fact that the phase has become trivial, one finds that it is still an interesting trivial phase, as was already observed in 18. Here, we study the zero modes of the model, and shed light on the zero mode present in one of the trivial phases.
By transforming c j → e ij ∇θ 2 c j as in the previous section, the Hamiltonian becomes As we show in the next section (where we consider the model for general parameters), the locations where the gap of this model closes are very similar to the locations of the phase transitions of the model with zero phase gradient, ∇θ = 0. Namely, they take the same form, if written in terms of the variablest = cos(∇θ/2)t and λ = cos(∇θ)λ, while µ remains unchanged. So, the gap closes whenλ = µ ±t, as well as when bothλ = −µ and |t| ≤ 2|µ|.
In Fig. 2, we show the phase diagram for a phase gradient ∇θ = π/3, and indicate the number of zero-modes on the left and the right separately. We first present the analysis of the zero modes and based on those results we explain the phase diagram of the model. Sticlet et al. 18 showed that the topological phase of this model has one zero mode at both edges as expected. The trivial phase, however, is divided into two regions. One region does not have any zero mode, while the other region has two 'Majorana' zero modes that are localized on one edge (i.e., a localized fermionic zero mode), while there is no zero mode on the other edge (see Fig. 2). The former trivial region corresponds to the trivial phase of the model without phase gradient while the later trivial region corresponds to the topological phase of the model without phase gradient with two Majorana zero modes on both sides. In what follows we present analytical wave functions for all the zero modes and determine for which parameters they are present. To find the Majorana zero modes we use Eqs. (14) and (15) and as before, we set g n = 1 2 (φ n + iψ n ). The 'bulk equations' read − µψ n + t cos( ∇θ 2 )ψ n−1 + λ cos(∇θ)ψ n−2 = 0, − µφ n + t cos( In this case, there are four boundary equations (two for each end) that differ from the bulk ones, namely and We start by solving the bulk equations, without paying attention to the boundary equations. We then solve the boundary equations, in the different regimes of the phase diagram.
The equations for φ n involves the solution for ψ n . Thus, the solution for φ n consists of two pieces, namely the general solution to Eq. (49) with the right hand side set to zero, which we will denote by φ gen,n , as well as a specific solution, for the full equation. We start with the ansatz φ gen,n ∼ x n . This gives us two correlation lengths 2λ cos(∇θ) .
(52) Thus, the generic solution is φ gen,n = L + x n + + L − x n − , where L ± are constants. As before,ψ n = φ gen,N −n+1 , which shows that the generic solution for φ is localized on the left edge and the solution forψ is localized on the right edge,ψ n = R + x N −n+1 , with R ± two constants. To find the full solution φ n , based on Eq.(49) we need to add a particular solution to φ gen,n of the form with constant S ± = κ ± R ± , since it should behave asψ. After some algebra, one finds that The general solution to the bulk equations (48) and (49) is thus given bỹ With the general solution for the bulk equations at hand, we turn our attention to the boundary equations, which we solve in the different regimes.
1) |x ± | > 1: In this case, both characteristic length scales are bigger than one, which occurs for the part of the phase diagram where the model without phase gradient is in the trivial phase. In this case, it is not hard to convince oneself that the boundary equations (50) and (50) lead to L ± = R ± = 0, which means that, as expected, there are no zero modes in this regime.
2) |x + | < 1 and |x − | > 1: In this case, the model is topological and x n − increases with n, which means that x n − is localized on the right edge instead of the left one. It is therefore convenient to write this solutions as this solution is localized on the right edge. The boundary equations (50) imply that R − = S − = 0. The boundary equations (51) give, after some algebra, thatL − = −R + t sin(∇θ/2) cos(∇θ)(µ+λ cos(∇θ)) , while S + = κ + R + as before. Thus, the solution for the zero mode is given bỹ We find that in this case, there is one zero mode, that is localized on both edges of the system. One special property of this zero mode, which differs from the case without a phase gradient, is that φ n has support on both edges of the system, while ψ n only has support on the right edge. Finally, we note that the case |x + | > 1 and |x − | < 1 is completely analogous.
3) |x ± | < 1: This case corresponds to the part of the phase diagram in which the model without phase gradient has two zero modes on both sides of the system. With the phase gradient, this model is in a trivial phase.
To determine if there are any zero modes, we again solve the boundary equations. The boundary equations for n = 1, 2, i.e. (50), give rise to terms that are proportional toψ at the left edge. Eq.(54) assures that these terms are of order x N ± and can be neglected in the thermodynamic limit. So the solution satisfies the boundary equations (50). The boundary equations (51) do give a non-trivial constraint. Namely, for a non-zero phase gradient ∇θ = 0 (for ∇θ = 0 the boundary equations are satisfied), one finds that These two boundary equations imply that R ± = 0. We conclude that in this regime there are two zero modes on the left side of the system, and none on the right side, i.e. φ n = L + x n + + L − x n − andψ n = 0. This precisely corresponds to the surprising result obtained by Sticlet et al 18 . We stress that this ordinary, or 'Dirac' zero mode on the left side of the system is not topological, but is in fact a consequence of fine tuning the parameters. We discuss this fine tuning in more detail in Sec. IV D. Nevertheless, as long as one keeps these parameters fine tuned, the only way this localized zero mode can disappear is via a closing of the gap, signaling a phase transition or a crossover.
We can now shed light on the phase diagram Fig. 2 (where we set ∇θ = π/3) and discuss it in more detail. The two solid lines in the figure, λ = µ/ cos(∇θ) ± t cos(∇θ/2)/ cos(∇θ), indicate phase transitions. Along the dashed line, λ = −µ/ cos(∇θ), the gap closes and a crossover occurs between a region with no zero mode and a region where the system has two zero modes on left edge and none on the right edge. The corners of the triangular region without any zero modes are given by (0, µ/ cos(∇θ)) and (±2µ/ cos(∇θ/2), −µ/ cos(∇θ)).
Upon increasing ∇θ, the slope of the lines indicating the phase transitions increases, and the size of the trivial center region increases. At ∇θ = π/2, the two phase transition lines are parametrized by t = ± √ 2µ and the trivial center region without any zero mode becomes a stripe in the middle of phase diagram. For ∇θ ∈ (π/2, π) the shape of the phase diagram is inverted with respect to the t-axis in comparison with Fig. 2.
We should note that the λ-axis (t = 0) should be treated separately. For t = 0, the model corresponds to two copies of the model that we studied in Sec. IV B, one for the even sites and another one for the odd sites. Therefore, the system is either in a trivial phase (for t = 0 and |µ| > |λ cos(∇θ)|) or in a topological phase (for t = 0 and |µ| < |λ cos(∇θ)|). In the topological phase the system has two Majorana zero modes on each edge, because the two chains are decoupled.
To provide further insight in the phase diagram (Fig. 2), we calculated the second derivative of the ground state energy for µ = 1, ∇θ = π/3 and t = 0 as well as t = 0.1 for λ ∈ [−4, 4] 2 . The result is shown in Fig. 3. We first consider the line t = 0 (the blue dashed line). As we described above, for λ = −4 the system is in the topological phase. The point λ = −2 is a critical point. Upon increasing λ, one enters the trivial phase. By passing the other critical point, namely λ = +2, one enters another topological phase. These two critical points give rise to divergencies in −d 2 E/dλ 2 as is shown in Fig. 3. This is a clear signature of a second order phase transition. For t = 0, (in the figure we used t = 0.1, the black line), the situation is quite different. At λ = −2, −d 2 E/dλ 2 is smooth (the critical point at λ = −2 and t = 0 shows its presence via a bump) and the system undergoes a crossover, although the gap closes. One still observes two divergencies symmetrically around λ = 2. These two divergencies correspond to the phase transitions indicated by the solid lines in Fig. 2.
We close this section by mentioning that it is of course possible to have the localized Dirac zero mode at the other edge of the system. One way to achieve this is by changing the phase dependence of the original pairing terms in the original Hamiltonian Eq. (46) to . This is equivalent to an inversion accompanied by the change ∇θ → −∇θ. We note that merely changing ∇θ → −∇θ does not change the role of the left and right hand side of the system. Basically the same calculation as above shows that this new pairing leads to two zero modes on right edge and none on the left edge (because the role of φ n and ψ n in Eqs. (48) and (49) is swapped). The model with these pairing terms has the same topological and trivial phases, however the left and right sides of the chain change their role. As we show in the next subsection, in the translational invariant formulation of the model, as in Eq. (47), the location of the zero-modes is determined by the relative sign of the phases of the hopping and paring terms.

D. The general case
To understand the fine tuning that is necessary to have the Dirac zero mode that resides on one side of the system as described in the previous section, we look at the more general Hamiltonian, where t 1 , ∆ 1 , t 2 and ∆ 2 are arbitrary complex parameters. In the case of periodic boundary conditions, we can define Ψ k = (c k , c † −k ) T and write the Hamiltonian as: where 1 is the two by two identity matrix. Performing a unitary transformation with U = 1 √ 2 (τ x + τ z ) we get, (60) By comparing the model we discuss here, Eq. (47) with (59) we find that in this case, all the imaginary parts depends are proportional to sin(∇θ/2) or sin(∇θ), for the nearest neighbor or NNN case respectively. So, these terms vanish for ∇θ = 0. In that case, we obtain Since we performed an unitary transformation, Det Q k = Det H k . So in the gapped phase, either topological or trivial, DetQ k = |h 3 (k) + ih 2 (k)| 2 = 0 can be used to define a topological invariant via the winding of Arg(h 3 (k) + ih 2 (k)), see 35 . This calculation leads to the same phase diagram we discussed before, see Fig. 1. We now consider a phase gradient, i.e. ∇θ = 0, which is the case we are interested in. Based on Eq. (47) we have The fact that (t 1 ) + (∆ 1 ) = (t 2 ) + (∆ 2 ) = 0, gives rise to Q k,11 = 0. This means that, similar to the ∇θ = 0 case, we have that DetQ k = |h 3 (k)+ih 2 (k)| 2 , despite the fact that Q k,22 = 0. Thus we find effectively the same phase diagram for the model with the phase gradient, namely the one given in Fig. 1, if we replace t → t cos( ∇θ 2 ) and λ → λ cos(∇θ). As we indicated in the previous section, by changing the paring terms in the original Hamiltonian Eq. (46) to , we can have the situation that the system has two 'Majorana' zero modes on the right edge and none on the left side. The gauge transformationc j = e ij ∇θ 2 c j changes these terms to jc j+1 + λ 2 e i∇θc † jc j+2 as before. In this case we find that (t 1 )− (∆ 1 ) = (t 2 )− (∆ 2 ) = 0, which results in Q k,22 = 0.
In class BDI, all the information about the zero modes is encoded in DetQ k . The discussion above shows that this is not so in the present case. Whether Q k,11 = 0 or Q k,22 = 0 plays an important role in determining the position of the (non-topological) localized zero modes. It is also clear what fine tuning we need in order to have a pair of 'Majorana' zero modes localized at one side of the system and none on the other. We need either Q k,11 or Q k,22 to be zero, but not both. This is the case if we fine tune (t 1 + ∆ 1 ) = (t 2 + ∆ 2 ) = 0 or (t 1 − ∆ 1 ) = (t 2 − ∆ 2 ) = 0, but not both which would imply that all these parameters are real, and one has an equal number of zero-modes on either side of the system.
To explore this situation further, we write the Hamiltonian in terms of Majorana operators, As a simple example we can see that for the Hamiltonian presented in Eq. (47), setting µ = 0 and ∇θ = π, yields following Majorana representation where it is evident that γ A,1 and γ A,2 do not appear in the Hamiltonian and therefore commute with it. Hence there are two Majorana zero modes on the left edge. It is interesting to note that for ∇θ = π, the Hamiltonian belongs to class BDI. From the form of the Hamiltonian in Eq. (47) this is not obvious, but it is for the form Eq. (46), because all the coupling constants are real. On the other hand, in this form, some of the couplings are staggered, and in the periodic case, the model is translationally invariant with a two-site unit cell. The phase diagram has a different structure in this case, with only three phases. The phase boundaries do not depend on t, and are given by λ = ±µ. Because of the two-site unit cell, we use the formulation of the phase-winding invariant as given by 39 . One finds that all three phases are in fact trivial. In the trivial phases with |λ| > |µ|, there is a localized Dirac zero mode only on the left side of the system, and no zero modes on the right side. This is consistent with the analysis of the model based on Eq. (47). From the point ∇θ = π it is clear that also in symmetry class BDI, there are Hamiltonians that have trivial phases, which have a localized Dirac zero mode only on one side of the system, if parameters are fine-tuned.
Our previous discussion led us to conclude that (t 1 + ∆ 1 ) = (t 2 + ∆ 2 ) = 0 could result in two zero modes on the left edge. Based on Eq. (63) we can see that this means that there should not be any terms like iγ A,j γ A,j+1 and iγ A,j γ A,j+2 present in the Hamiltonian. We can shed more light on this issue based on our analytical solution for the non-uniform pairing with nearest neighbor hopping and pairing.
As a first a step, we assume that t = 0. This means that we have two decoupled chains with a phase gradient. Our previous analysis shows that in the topological phase we have one Majorana zero mode on each edge. The wavefunctions for these Majorana modes are given in Eqs. (44) and (45). The crucial difference between these two wave functions originates in the direction of the phase gradient, which causes the left mode to be independent of γ B , while the right modes consists of both γ A and γ B . Namely, for the left mode g n is purely real while for the right mode g n is complex, and hence involves both φ andψ.
In the second step, we turn on nearest neighbor couplings, i.e. t = 0. We see that the first four terms in the Hamiltonian Eq. (63) result in a coupling between the zero modes of the two decoupled chains. Under the assumption that (t 1 + ∆ 1 ) = 0 and (t 2 + ∆ 2 ) = 0, which holds in our analytic calculation of the zero modes, we find that the right zero modes from the two different chains become coupled because of the iγ B,j γ B,j+1 terms present then t = 0, which gaps them out. On the other hand, the zero modes on the left edge do not become coupled directly, and remain gapless. Their wavefunctions are modified to the new ones presented in Eq. (54).
Finally, we mention that we checked numerically that under the conditions (t 1 + ∆ 1 ) = (t 2 + ∆ 2 ) = 0, the system has a phase for which the ground state has two zero modes located on the left edge, and none on the right edge. The same holds true in the case that (t 1 − ∆ 1 ) = (t 2 − ∆ 2 ) = 0, if one exchanges the left and right edge of the system.

V. DISCUSSION
In this paper, we investigated the 'one-sided' fermionic zero modes observed by Sticlet et al. 18 , by solving the Kitaev model, in the presence of complex hopping and pairing terms, including NNN terms, for open chains. From our investigation, it became clear that fine-tuned parameters are necessary for such zero modes to exist, but under the fine-tuned conditions, the gap needs to close in order to destroy them. Leaving the fine-tuned conditions gaps these zero-modes out, turning them into one-sided low-energy subgap modes. Phases with such one-sided bound states can occur both in one-dimensional systems in class D, as well as in class BDI. These modes are not protected by topology, which means that they can occur in the topologically trivial phase.
The general condition for the existence of 'one-sided zero modes' is most easily explained in terms of the Majorana formulation of the chains. Starting from a situation in which two pairs of (delocalized) Majorana bound states are present (i.e., in class BDI), one needs a perturbing term such that the two Majoranas describing the mode on, say, the left side are coupled, while the modes on the right side are not. It is worth to mention that if one assumes a different phase gradient for the nearest neighbour and next-nearest neighbour pairing terms in Eq. 46, these 'one-sided zero modes' gap out.
There has been a lot of progress on models in higher dimensions, that exhibit exact zero modes, see for instance Ref. 40. It would be interesting to investigate if it is possible to construct models, that exhibit 'one-sided' zero modes along the lines of the ones described in this paper, even in those higher-dimensional systems.
Acknowledgements -We would like to thank F. Pollmann, C. Spånslätt and R. Verresen for interesting discussions. This research was sponsored, in part, by the Swedish research council. In this appendix we present details of the solution for the open Kitaev chain with generic real parameters and free boundary conditions. The Hamiltonian reads: It is helpful to recall the (of course well known) solution for the periodic case, which is obtained via a Fourier transformation, c j = 1 √ N k e ikj c k , and defining Ψ k = (c k , c † −k ) T . This results in Diagonalization of this 2 × 2 matrix gives us: where f k is a new fermionic quasiparticle annihilation operator.
To tackle the open case, we use the LSM method which is reviewed in Sec. II. To this end, we need to arrange the Hamiltonian to have the form of Eq. (1), To find φ α and ψ α from Eqs. (10) and (11) i.e., α ψ α , we have to construct the matrices A − B and A + B.
We present these matrices for the more general case of Hamiltonian in Eq. (58), i.e., because we need them later on. In this case, A − B and A + B read, For the Hamiltonian Eq. (A1), i.e. with t 1 = 1, ∆ 1 = ∆, and t 2 = ∆ 2 = 0, these reduce to, Using these matrices in Eq. (10) one gets for 3 ≤ n ≤ N − 2. We call this the 'bulk equation'. In the case of periodic boundary conditions, this is actually the only equation one has to consider. However, for an open chain with free boundary conditions, we also have four boundary equations which are different from the bulk one, namely for n = 1, 2, N − 1 and N one has: We note the difference between the φ α,1 term in the equation for n = 1 and the φ α,N term in the equation for n = N .
To solve these equations we can start with an ansatz for the eigenvalues Λ α . Note that the bulk equation is the same for both the periodic and the open chain. This suggests to use our knowledge about the periodic case.
Now we need to find an equation based on which one can determine all the possible values of α. With the ansatz for Λ α , we solve Eq. (A11) by the standard approach, i.e. we consider φ α,n ∼ x n α . Using this in Eq. (A11) gives us: One checks that e ±iα are solutions independent of the parameter K. Since we have found two roots, we can find the other two, which are given by e ±iβ where β satisfies cos α + cos β = K 2 . (A19) Therefore each α has a β partner. We note that α and β are equivalent. The associated eigenvalues can be written in the same functional form, i.e. Λ α = Λ β = (µ − cos β) 2 + ∆ 2 sin 2 β, which follows from Eq. (A19). We continue to use α as the label indicating the eigenvalue.
These solutions tell that e ±inα and e ±inβ are the most general solution for the bulk equation. Now we need to determine a linear combination of these functions that satisfies the boundary equations. Treating the left and right edges in an equivalent way, we consider the following combination: in which A 1 ,A 2 ,B 1 and B 2 are constants. Using this ansatz, Eqs. (A13) and (A14) give us: Based on these relations, we rewrite the ansatz: Finally, we make sure that the ansatz satisfies Eqs. (A12) and (A15), which leads to the following equations: in terms of the functions To find a non-trivial solution for A 1 and A 2 , we require that the determinant of the matrix in Eq. (A24) is zero. This gives us another equation for α and β: This equation should be solved together with Eq. (A19) to give us all admissible labels. Generically, this has to be done numerically.
In the analysis below, we focus on the regime with µ ≥ 0 and ∆ ≥ 0. We assume that ∆ = 1, the case ∆ = t = 1 was considered explicitly in 20,21 . From the equations (A28) and (A19), we see that a solution (α, β) for ∆ > 0 also gives a solution for ∆ < 0 (though the form of the wave function φ α,n changes). In addition, the solutions for µ < 0 can be related to the solutions with µ > 0. If a pair (α, β) satisfies the equations for µ > 0, the pair (α + π, β + π) will satisfy the equations for µ < 0. Note that this shift does not change Eq. (A28). However, it gives rise to a minus sign in the left hand side of Eq. (A19) which indeed changes the sign of µ. Finally, the actual eigenvalues Λ α are also unchanged.
Thus from now on, we assume that µ, ∆ ≥ 0. The structure of the solutions (α, β) is as follows. For µ > 1, one finds N solutions, for which α and or β is real. Because α and β are completely equivalent, we assume that α is real. When 0 ≤ µ < 1, there are N −1 solutions, with α real, and β either real or complex. We note that if β is complex, its real part Reβ = 0 for ∆ < 1, and Reβ = π for ∆ > 1. The 'missing' solution has both α and β complex, and corresponds to the zero mode, which we describe in detail below. In Fig. 4, we show this for a chain of N = 6 sites, ∆ = 0.8 and different values of µ.
Before we do so, we first discuss the solutions with α real. We first note that for any (α, β) pair that solves Eqn. (A28) and (A19), all the combinations of (±α, ±β) are also a solution. Since these pairs give rise to same wavefunction, we only consider α in the range 0 ≤ α ≤ π.
The solutions are then obtained by finding the solutions of Eq. (A28), where β is given by Eq. (A19). Special care has to be taken in the case that both α and β are real, say (α, β) = (α 1 , β 1 ), because one will also find the equivalent solution (α, β) = (β 1 , α 1 ), so one has to restrict the range of α further, to avoid 'double counting' of solutions.
Nevertheless, the value α c is useful when specifying the appropriate range for α. If there are solutions with both α and β real, one has that either α < α c < β, or β < α c < α. In addition, for the range ∆ ≤ √ 1 − µ, one finds that all the solutions (α, β) with β imaginary have α > α c . Thus, to find all solutions in this range, one should only take the solutions for α such that α c < α < π. For the range ∆ ≥ √ 1 + µ, the situation is opposite, and one should take the solutions for α in the range 0 ≤ α < α c . In the other regime, namely √ 1 − µ < ∆ < √ 1 + µ, one has to consider all solutions for α in the range 0 ≤ α < π.
We now turn our attention to the Majorana zero mode solution. The goal is to find the analytical expression for the wave function of this mode. For simplicity, we work in the limit of large system size, i.e., N → ∞.
By analyzing Eq. (A28), one finds that the solution one loses, is the one with smallest positive, real α. Taking the limit α → 0 and N → ∞ of Eq. (A28), using Eq.(A19), gives This shows that there is a solution with α = 0, for µ = 1.
In addition, further analysis shows that for µ < 1, one loses this solution, both for ∆ < 1 and ∆ > 1, while for µ > 1, this solution shifts to finite, positive α. This behavior can be seen for a chain with N = 6 sites, ∆ = 0.8 and µ = 1.2, 0.6, 0.25 in Fig. 4. In the case of µ = 0.25, only the solutions with α > α c ≈ 0.25π are independent, so the there are still only five solutions. The additional, sixth solution is still a zero-mode. We note that for finite N , the value of µ for which one loses the solution has 1/N corrections, and depends on ∆. That the phase transition between the trivial and topological phase occurs for µ = 1 in the large N limit is of course well known, and is given by the value of µ for which the gap closes. Based on Eq. (A4), we infer that µ = ±1 are the only possible values of chemical potential for which gap closes (provided that ∆ = 0). Now we turn to finding the missing root and its associated features. To do so we need to consider different cases.
1) ∆ < 1 and √ 1 − ∆ 2 < µ < 1: In this regime, we lost one solution with α real, so we look for a solution with both α and β imaginary, and in fact, purely real. Such a solution indeed exist namely, which solves Eq. (A19) and Eq. (A28) in the large N limit. For √ 1 − ∆ 2 < µ < 1, both ξ 1 and ξ 2 are real. Let us explore the properties of this solution. First, substituting this result back into the Eq. (A16) gives us Λ α * = 0, so we indeed have a zero-mode. This means that we can use Eq. (8) to solve for the wave function. Alternatively, we can set A 1 = 0 in Eq. (A23) to obtain the Majorana mode that is localized on the left side of the system. Either approach gives φ α * ,n = Ce − n ξ 1 sinh( n ξ 2 ), where C is a normalization constant. Because ξ 1 < ξ 2 , the mode φ α * is indeed localized on the left edge. The same reasoning can be done for ψ α,n . The important observation is that (A + B)(A − B) has the same structure as (A − B)(A + B) if we look at it from the other side of the chain, i.e. n → N + 1 − n. So we get ψ α * ,n = φ α * ,N +1−n , which tells us that ψ α * ,n is localized on the right edge.
2) ∆ < 1 and µ < √ 1 − ∆ 2 : For µ < √ 1 − ∆ 2 , the parameter ξ 2 in Eq. (A30) becomes imaginary, so is more natural to rewrite the previous solution. The root can be written as: Again, one finds that Λ α * = 0. Using the same logic as above, one finds that φ α * ,n = Ce − n ξ sin(nq), with C some constant. This result shows that φ α * is localized on the left edge. Although this in this case instead of having decaying functions, we have an oscillatory decaying function.
3) ∆ > 1: In this case we can not use the previous results, because √ 1 − ∆ 2 becomes imaginary. One finds that the new root in this regime is given by We see that ξ 1 < ξ 2 since µ < 1. One can check that for this root Λ α * = 0, hence it is also a zero mode. To find the Majorana mode that is localized on the left edge, we again set A 1 = 0 in Eq. (A23), which results in φ α * ,n = Ce This result shows that φ α * is localized on the left edge.
To close this section we note that for µ = √ 1 − ∆ 2 , we have α * = β * . Therefore one can not use x n α and x n β as separate solutions, but one should use nx n α as the other independent solution.
Appendix B: The zero-modes of the Kitaev chain with a phase gradient In this appendix, we investigate the zero mode of the Kitaev chain, in the presence of a phase gradient in the order parameter. We assume that |t| = |∆| = 1.
As we mentioned in Sec. IV B, after a gauge transformation the Hamiltonian takes the form in which ∇θ is the phase gradient per site, which is constant. To find the zero-mode, we use the method which is presented in Sec. II. From Eq. (A7), the matrices A − B and A + B read (B3)