Exact eigenvectors and eigenvalues of the finite Kitaev chain and its topological properties

We present a comprehensive, analytical treatment of the finite Kitaev chain for arbitrary chemical potential and chain length. By means of an exact analytical diagonalization in the real space, we derive the momentum quantization conditions and present exact analytical formulas for the resulting energy spectrum and eigenstate wave functions, encompassing boundary and bulk states. In accordance with an analysis based on the winding number topological invariant, and as expected from the bulk-edge correspondence, the boundary states are topological in nature. They can have zero, exponentially small or even finite energy. Further, for a fixed value of the chemical potential, their properties are ruled by the ratio of the decay length to the chain length. A numerical analysis confirms the robustness of the topological states against disorder.


Introduction
The quest for topological quantum computation has drawn a lot of attention to Majorana zero energy modes (MZM), quasiparticles obeying non-Abelian statistics hosted by topological superconductors [1]. The archetypal model of a topological superconductor in one dimension was proposed by Kitaev [2]. It consists of a chain of spinless electrons with nearest neighbour superconducting pairing, a prototype for p-wave superconductivity. As shown by Kitaev in the limit of an infinite chain, for a specific choice of parameters, the superconductor enters a topological phase where the chain can host a couple of unpaired zero energy Majorana modes at the end of the chain [2]. This model has thus become very popular due to its 3 Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. apparent simplicity and it is often used to introduce topological superconductivity in one dimension [1,3]. Also more sophisticated realizations of effective p-wave superconductors, based on semiconducting nanowire-superconductor nanostructures [4][5][6][7][8][9][10][11], ferromagnetic chains on superconductors [12][13][14][15][16] or s-wave proximitized carbon nanotubes [17][18][19], all rely on these fundamental predictions of the Kitaev model. While the theoretical models are usually solved in analytic form for infinite or semi-infinite chains, the experiments are naturally done on finite-length systems. For example, for the iron chain on a superconductor investigated in [13] it is expected that the chain length is shorter than the superconducting coherence length [20]. Spectral properties of a finite-length Kitaev chain have been addressed in more recent papers [21][22][23][24], and have confirmed the presence of bound states of exponentially small energy in sufficiently long finite Kitaev chains.
As noticed by Kao using a chiral decomposition [21], a finite-length Kitaev chain also supports modes with exact zeroenergy. However, they are only found for discrete values of the chemical potential. As shown by Hegde et al, in [22], these exact zeros can be associated to a fermionic parity crossing in the open Kitaev chain. Investigations of the finite chain have also been performed by Zvyagin [23] using the mapping of a Kitaev chain onto an X-Y model for N spin 1/2 particles in transverse magnetic field, for which convenient diagonalization procedures are known [25,26]. Kawabata et al in [27] demonstrated that exact zero modes persist also in a Kitaev chain with twisted boundary conditions. Despite much being known by now about the finite length Kitaev chain, and in particular its low energy properties, the intricacy of the eigenvector equation still constitutes a challenge. In this work we address this longstanding problem. Exact results for the full energy spectrum and the associated bound states are provided for arbitrary chemical potential by analytical diagonalization in the real space. Our results are not restricted to long chains or to long wavelengths, and thus advance some of the findings in [22,23] and respectively [24]. They also complete the analysis of the eigenstates of an open Kitaev chain performed by Kawabata et al [27] which was restricted to the exact zero modes. This knowledge allows one a deeper understanding of the topological properties of the Kitaev chain and, we believe, also of other one-dimensional p-wave superconductors.
Before summarizing our results, we clarify the notions used further in our paper. (i) Since 'phase' properly applies only to systems in the thermodynamic limit, we shall use 'topological regime' to denote the set of parameters for which a topological phase would develop in an infinite system. (ii) We shall call 'topological states' all boundary states of the finite system whose existence in the topological regime is enforced by the bulk-boundary correspondence [1]. Hence, the existence of the topological states is associated to a topological invariant of the bulk system being non trivial. Quite generally, the topological states can have zero, exponentially small, or even finite energy. (iii) When the gap between the topological states and the higher energy states is of the same order or larger than the superconducting gap Δ, we consider them to be 'robust' or 'topologically protected'. Their energy may be affected by perturbations, but not sufficiently to make them hybridize with the extended (bulk) states. (iv) All Hamiltonian eigenstates can be written as superpositions of Majorana (selfconjugate) components. When the energy of the topological state is strictly zero, the whole eigenstate has the Majorana nature and becomes a true 'Majorana zero mode' (MZM).
Using our analytical expressions for the eigenstates of a finite Kitaev chain, we recover in the limit of an infinite chain the region for the existence of the MZM given by the bulk topological phase diagram; the latter can be obtained using the Pfaffian [2] or the winding number [28] topological invariant. For a finite-length chain MZM only exist for a set of discrete values of the chain parameters, see equation (87) below, in line with [21,27]. These states come in pairs and, depending on the decay length, they can be localized each at one end of the chain but they can also be fully delocalized over the entire chain. Even in the latter case the states are orthogonal and do not 'hybridize' since they live in two distinct Majorana sublattices. Similar protection of topological zero energy modes living on different sublattices has recently been observed experimentally in molecular Kagome lattices [29].
The paper is organized as follows. Section 2 shortly reviews the model and its bulk properties. Section 3 covers the finite size effects on the energy spectrum and on the quantization of the wave vectors for some special cases, including the one of zero chemical potential. The spectral properties at zero chemical potential are fully understood in terms of those of two independent Su-Schrieffer-Heeger-like (SSHlike) chains [27,30]. The eigenstates at zero chemical potential, the symmetries of the Kitaev chain in real space, as well as the Majorana character of the bound state wave functions are discussed in section 4. In sections 5-7 we turn to the general case of finite chemical potential which couples the two SSH-like chains. While section 5 deals with the energy eigenvalues and eigenvectors of the finite chain, section 6 provides exact analytical results for the MZM. In section 7 the influence of disorder on the energy of the lowest lying states is investigated numerically. In section 8 conclusions are drawn. Finally, appendices A-G contain details of the factorisation of the characteristic polynomial in real space and the calculation of the associated eigenstates.

Model
The Kitaev chain is a one dimensional model based on a lattice of N spinless fermions. It is characterized by three parameters: the chemical potential μ, the hopping amplitude t, and the p-wave superconducting pairing constant Δ. The Kitaev Hamiltonian, written in a set of standard fermionic operators d j , d † j , is [1,2] (1) where the p-wave character allows interactions between particles of the same spin. The spin is thus not explicitly included in the following. We consider Δ and t to be real parameters from now on.
The Hamiltonian in equation (1) has drawn particular attention in the context of topological superconductivity, due to the possibility of hosting MZM at its end in a particular parameter range [2]. This can be seen by expressing the Kitaev Hamiltonian in terms of so called Majorana operators γ A,B , yielding the form Notice that, in virtue of equation (2) it holds {γ A j , (γ A j ) † } = 2 (γ A j ) 2 = 1, and similarly for γ B j . For the particular parameter settings Δ = ±t and μ = 0, which we call the Kitaev points, equation (3) leads to a 'missing' fermionic quasiparticle q ± : This quasiparticle has zero energy and is composed of two isolated Majorana states localised at the ends of the chain. In general, the condition of hosting MZM does not restrict to the Kitaev points (μ = 0, Δ = ±t). Further information on the existence of boundary modes is evinced from the bulk spectrum and the associated topological phase diagram. If the boundary modes have exactly zero energy, their Majorana nature can be proven by showing that they are eigenstates of the particle-hole operator P. Equivalently, if γ † M is an operator associated to such an MZM it satisfies γ † M = γ M and γ 2 M = 1/2. The topological phase diagram is shortly reviewed in section 2.3.

Bulk spectrum
The Hamiltonian from equation (1) in the limit of N → ∞ reads in k spacê where we introduced the operators d k = 1 √ N j e −i j kd d j and k lies inside the first Brillouin zone, i.e. k ∈ − π d , π d and d is the lattice constant. The 2 × 2 Bogoliubov-de Gennes (BdG) matrix is easily diagonalized thus yielding the excitation spectrum Note that for μ = 0 equation (7) predicts a gapped spectrum whose gap width is either 4Δ (|Δ| < |t|) or 4t (|t| < |Δ|).

Topological phase diagram
The BdG Hamiltonian (6) is highly symmetric. By construction it anticommutes with the particle-hole symmetry P = σ x K, where K accounts for complex conjugation. The particle-hole symmetry turns an eigenstate in the k space corresponding to an energy E and wavevector k into one associated with −E and −k. The time reversal symmetry is also present in the Kitaev chain and is given by T = K. Finally, the product of T P = C = σ x is the chiral symmetry, whose presence allows us to define the topological invariant in terms of the winding number [28]. Note that all symmetries square to +1, placing the Kitaev chain in the BDI class [31]. The winding number is given by [28,32] where w(k) = arg [2Δ sin(kd) + i (μ + 2t cos(kd))] and ∂ k w(k) is the winding number density. A trivial phase corresponds to ν = 0, a non trivial one to finite integer values of ν. The winding number relates bulk properties to the existence of boundary (not necessarily MZM) states in a finite chain. This property is known as bulk-edge correspondence. Due to their topological nature, their existence is robust against small perturbations, like disorder. This point is further discussed in section 7. The phase diagram constructed using the winding number invariant is shown in figure 1. The meaning of two different values for the winding number is clearer when we recall the Kitaev Hamiltonian in the Majorana basis. In a finite chain the leftmost lattice site consists of the A Majorana operator γ A 1 connected to the bulk by the i(Δ − t) hopping and the B Majorana operator γ B 1 connected by the i(Δ + t) hopping. With Δ > 0 and t < 0 (the ν = +1 phase) the Majorana state at the left end of the chain will consist mostly of the weakly connected γ B 1 . If t > 0 (the ν = −1 phase), γ A 1 is connected to the bulk more weakly and contributes most to the left end bound state.
The boundaries between different topological phases can be obtained from the condition of closing the bulk gap, i.e. E ± (k) = 0 (cf equation (7)). That is only possible if both terms under the square root vanish. The condition of Δ = 0 forces the gap closing to occur at kd = 0 or k = πd, and the remaining term vanishes at these momenta if μ = ±2t. The four insets in figure 1 show the behavior of w(k), leading to either a zero (for w(−π) = w(π)) or non zero winding number, see equation (8).
Physically speaking, the Kitaev chain is in the topological phase provided that Δ = 0 and the chemical potential lies inside the 'normal' band (|μ| 2|t|).

Spectral analysis of the finite Kitaev chain
One of the characteristics of finite systems is the possibility to host edge states at their ends. To account for the presence and the nature of such edge states, we consider a finite Kitaev chain with N sites and open boundary conditions, yielding N allowed k values. In this section we shall consider the situation in which one of the three parameters Δ, t and μ is zero. Already for the simple case μ = 0 and Δ = 0, t = 0 the quantization of the momentum turns out to be non trivial. The general case in which all parameters are finite is considered in sections 5-7.
We start with the BdG Hamiltonian of the open Kitaev chain in real space, expressed in the basis of standard fermionic Then where the BdG Hamiltonian H KC is These matrices have the tridiagonal structure The spectrum can be obtained by diagonalisation of H KC in real space. We consider different situations.

Δ = 0
The BdG Hamiltonian is block diagonal and its characteristic polynomial P λ (H KC ) = det [λ − H KC ] factorises as The tridiagonal structure of C straightforwardly yields the spectrum of a normal conducting, linear chain [33] E Δ=0 where n runs from 1 to N. Since k n ∈ R, only bulk states exist for Δ = 0.

t = 0
In the beginning we consider both t and μ to be zero and include μ = 0 in a second step. The parameter setting leads to a vanishing matrix C, see equation (11), and the characteristic polynomial of the system reads: where we used the property S † = −S. Due to the fact that the commutator [ , S] = 0 vanishes 4 , one finds [34] The characteristic polynomial can still be simplified to the product The matrix iS is hermitian and describes a linear chain with hopping iΔ. As a consequence, we find the spectrum to be [33] where n runs from 1 to N and each eigenvalue is twice degenerated. Notice the phase shift by π/2 compared to the spectrum of an infinite chain equation (7). We discuss this phase shift in more detail in section 3.3. Furthermore if, and only if, N is odd, we find two zero energy modes, namely for n = (N + 1)/2. Their existence for odd N and the degeneracy is due to the chiral symmetry.
The chemical potential μ can be included easily. Exploiting the properties of H KC , we find the characteristic polynomial to be with Λ 2 := λ 2 − μ 2 . The same treatment as in the previous where n runs again from 1 to N. Again no boundary modes are found for t = 0. 4 Note that λ can be zero and it will for odd N. Hence, the standard formula to calculate the determinant of a partitioned 2 × 2 matrix can not be used here, because it requires the inverse of one diagonal block. We use instead Silvester's formula [34]: det

μ = 0
The calculation of the spectrum for μ = 0 requires a more technical approach, since the structure of the BdG Hamiltonian equation (10) prohibits standard methods. One important feature of the Kitaev chain can be appreciated inspecting equation (3). The entire model is equivalent to two coupled SSH-like chains [30,35] containing both the hopping parameters a := i (Δ − t) and b := i (Δ + t), see figure 2. Explicitly, where N 1,2 depend on N. If N is even we have N 1 = N/2 and Independent of the number of atoms, the first and the second lines in equation (21) describe two SSH-like chains, coupled by the chemical potential μ. We define here the SSH-like basis of the Kitaev chain as: where '|' marks the boundary between both chains. We call the first one, starting always with γ A 1 , the α chain, and the second the β chain, such that Ψ even,odd Hamiltonian in the SSH-like basis reads The independent SSH-like chains are represented by the square matrices H α and H β of size N. Both chains are coupled by the matrices τ and τ † , which contain only the chemical potential μ, in a diagonal arrangement specified below.
The pattern of these matrices is slightly different for even and odd number of sites. If N is even we find H even and τ even = −iμ N/2 ⊗ τ z , where τ z denotes the Pauli matrix. The odd N expressions are achieved by removing the last line and column in H even α , H even β and τ even . As shown in more detail in appendix B, for μ = 0 the characteristic polynomial can be expressed as the product of two polynomials of order N where the product form reflects the fact that the Kitaev chain is given in terms of two uncoupled SSH-like chains, as illustrated in figure 2. Even though the polynomials ζ N and N belong to different SSH-like chains, both obey a common recursion formula typical of Fibonacci polynomials [36][37][38] and differ only in their initial values ⎛ Fundamental properties of Fibonacci polynomials are summarized in appendix A. The common sublattice structure of both chains sets the stage for a relationship between ζ j and j : the exchange of a's and b's enables us to pass from one to the other Moreover, equation (27) implies that a Kitaev chain with even number of sites N is fundamentally different from the one with an odd number of sites. This property is a known feature of SSH chains [39]. The difference emerges since, according to equations (B20), (B22), it holds because the number of a and b type bondings in both subchains is the same. This leads to twice degenerate eigenvalues. An equivalent relationship for even N does not exist. The closed form for ζ j and j , as well as their factorization, is derived in appendix B. The characteristic polynomial can be used to obtain the determinant of the Kitaev chain, here for μ = 0, because evaluating it at λ = 0 leads to: According to equation (26) we need only to know ζ N and N at λ = 0. The closed form expression for ζ j at λ = 0 reduces to while j | λ=0 . follows from equation (29). We find that there are always zero energy eigenvalues for odd N, but not in general for even N, as it follows from det H μ=0 Additional features of the spectrum are discussed in the following.
3.3.1. Odd N. The spectrum for odd N is given by two contributions E μ=0 where k n d = nπ/(N + 1) and n runs from 1 to N, except for n = (N + 1)/2. This constraint on n is a consequence of the equations (67), (69) below which show that the boundary condition equation (70) cannot be satisfied for kd = π/2. Hence, no standing wave can be formed. Each zero eigenvalue belongs to one chain. As discussed below, two decaying states are associated to equation (33), whose wave functions are discussed in section 4.2. These states are MZM.
The zeros of f β,α (k) define the proper wave vectors k 3,2,1 of the finite system and these cut the dispersion relation at the correct positions, such that E ± (k j ) = ±E j . (b) Zoom of (a) for k ∈ 0.5/d, where the momenta k are in general not equidistant in the first Brillouin zone. Rather, the bulk quantization condition follows from the interplay between Δ and t and is captured in form of the functions f β,α (k) (cf appendix B.2), whose zeros f β,α (k) define the allowed values of k. Note that kd = 0, π/2 are excluded as solutions, due to their trivial character. The functions f β,α (k) follow from the factorisation of the polynomials N and ζ N . The negative sign in equation (36) belongs to the α subchain, while the positive one to the β subchain. The spectrum following from equation (37) is illustrated in figure 3. We observe that equations (35) and (37) hold for all values of t and Δ, independent of whether |Δ| is larger or smaller than |t|. The two situations are connected by a phase shift of the momentum kd → kd + π/2, which influences both the spectrum and the quantization condition. In the end all different ratios of Δ and t are captured by equations (35) and (37), due to the periodicity of the spectrum. However, when we consider decaying or edge states this periodicity is lost (see equations (40) and (41) below) and |t| ≶ |Δ| lead to different quantization rules. The hermiticity of the Hamiltonian allows a pure imaginary momentum for μ = 0, but a simple exchange of k to iq in equation (36) does not lead to the correct results. We introduce here the functions h β,α (q) := tanh [qd (N + 1)] ± m tanh (qd) , (38) similar to f β,α (k) in equation (36), where m contains both ratios of Δ and t: Again, the positive sign in equation (38) belongs to the β chain and the negative one to the α chain. The exact quantization criterion is provided by the zeros of h β,α (q), as illustrated in figure 4. The associated energies follow from the dispersion relation We notice that equation (41) is only well defined for zero or positive arguments of the square root. Indeed, all solutions of equation (40), if existent, lie always inside this range, because using h β,α (q) = 0 in equation (41) yields Hence, each wavevector from equation (40) corresponds to two gap modes, since the gap width is 4 min {|Δ|, |t|} and the fraction inside equation (42) is always smaller than one. We can restrict ourselves to find only positive solutions qd, due to the time reversal symmetry. The number of physically different solutions of equation (40) is zero or two and it follows always from the equation containing the positive factor m or −m. Consequently, according to equation (38), only none or two gap modes can form and both belong to the same subchain, α or β. Moreover a solution exists if, and only if, |m| ∈ [1, N + 1].
In the limiting case when |m| → 1, i.e. at the Kitaev points, the solution qd → ∞ and the associated energies E ± from equation (42) go to zero. The eigenstate will be a Majorana zero energy mode, see section 4.1.
In the second special case of |m| → N + 1 the solution approaches zero. The value q = 0 is only in this particular scenario a proper momentum, see appendix B.2. The momentum q = 0 yields the energies E ± (0) = ±2 min {|Δ|, |t|}, which mark exactly the gap boundaries.
Increasing the value of |m| beyond N + 1 entails the absence of imaginary solutions. The number of eigenvalues of a Kitaev chain is still 2N for a fixed number of sites and consequently equation (37) leads now to N real values for kd, instead of N − 1. In other words, the two former gap modes have moved to two extended states and their energy lies now within the bulk region of the spectrum, even though the system is still fully gaped. This effect holds for the Kitaev chain as well as for SSH chains. Physically this means, that a 'boundary' mode with imaginary momentum q and corresponding decay length ξ ∝ 1/q reached the highest possible delocalisation in the chain.
The limit of N → ∞ yields always two zero energy boundary modes; since the momentum is qd = arctanh(|1/m|), due to equations (38), (40) and according to equation (42) the energy goes to zero. If we consider the odd N situation in the limit of an infinite number of sites, we have there two zero energy boundary modes as well. The results of this section are summarized in table 1.

Eigenvectors (μ = 0)
We use the SSH-like basis to calculate the eigenvectors of the Hamiltonian equation (23) at μ = 0. The eigenvectors ψ are defined with respect to the SSH-like chains α and β, see equation (23), with the feature that always either v β or v α can be chosen to be zero, yielding the solutions ψ α and ψ β , respectively We are left to find the eigenvectors of a single tridiagonal matrix which we did basing on, and extending the results of [40]. We focus here on the edge and decaying states, while the rest of our results is in appendix C. Remember that in the SSH-like basis equation (22) the Majorana operators γ A j and Table 1. Overview of the quantization rule for the wave vectors of the finite Kitaev chain in different scenarios. The wavevectors k, k n , κ 1,2 used together with equation (7) and q with equation (41) yield the correct finite system energies and k, k n , q ∈ R, κ 1,2 ∈ C.

Equation for eigenstate Majorana Requirements
Quantisation rule Zero modes elements character Yes, if for some n: j , alternate at each site, thus defining two interpenetrating 'A' and 'B' type sublattices.

Even N
We define the vectors v α and v β via the entries where x, y and X , Y are associated to the A and B sublattices, respectively. The internal structure of v α ( v β ) reflects the unit cell of an SSH-like chain and thus simplifies the calculation.
In the real space x l (X l ) belongs to site j = 2l − 1 and y l (Y l ) to j = 2l, where j = 1, . . . , N. Searching for solutions on the subchain α implies setting v β = 0 and solving H even and where l runs from 1 to N/2 − 1. The solution for Δ = ±t is (in agreement with [40]) x l where l = 1, . . . , N/2, θ/(2d) denotes the momentum k (q) for extended (gap) states and E ± is the dispersion relation associated to k (equation (35)), or q (equation (41)). The entries of the eigenvectors are essentially sine functions for the extended states and hyperbolic sine functions for the decaying states where the prefactor s depends on the ratio of Δ and t: An illustration of ψ α is given in figure 5. The allowed momenta k or q follow from the open boundary conditions The first condition is satisfied due to T 0 (θ) = 0 for any momentum. The second condition yields the quantization rules f α (k) = 0 and h α (q) = 0 for the α chain, see equations (37) and (40).
The eigenvector ψ β entails v α = 0 and the entries of v β follow essentially by replacing a's and b's in the equations (51) and (52). We find where j = 1, . . . , N/2 and Δ = ±t. The quantisation condition follows from the open boundary condition: . Further, from the quantization rules it follows that gap modes belong always to the same subchain α or β for even N.  (57) , which holds for all eigenstates. Together with equations (51), we find in general and similarly X N Recalling the definition of the SSH-like basis, equation (22), and introducing the operators ψ † α,β associated to the states ψ α,β in equation (44), we find the expression where v α is the norm of the vector v α . A similar term is found for ψ † β . We notice that equations (59) and (60) below are true for all kinds of eigenstates, i.e. extended, decaying states and MZM, of the BdG Hamiltonian in equation (23) at μ = 0. The character (statistics) of the operators depends on The symmetry in equations (58) states that (ψ † α ) 2 is essentially determined by (x N/2 /y 1 ) 2 . For a = 0 and consequently E = 0, we find from equations (47), (48) and (58) that (x N/2 /y 1 ) 2 = Notice that for the chosen parameter set the gap state is more localized than the one in figure 5. In contrast the extended state has lower weight at the ends of the chain. The gap mode (bulk state) is associated with q = 0.2027/d (k = 1.4886/d) and Thus the operators associated to the finite energy states ψ α , including the ones depicted in figures 5(a) and 6(a), obey fermionic statistics. This result holds also true in the case a = 0 and E = 0 as can be seen by using the corresponding eigenstates (appendix C). Similar results hold for (ψ † β ) 2 . We turn now to Majorana zero modes, which at μ = 0 only exist at the Kitaev points Δ = ±t, see equation (32).
When Δ = t we find two zero energy modes ψ α and (ψ † α ) 2 = 1/2 in equation (60). In contrast, both zero energy modes are on the β chain for Δ = −t. We find ψ β These states are the archetypal Majorana zero modes [1,2]. Due to their degeneracy, these modes can be recombined into fermionic quasiparticles by appropriate linear combination, see in equations (4a) and (4b) from section 2.1.
The difference to the edge states on subchain α is the exchanged role of Majorana operators γ A j , γ B j . The chosen parameters are t = −5 meV and Δ = 1 meV. The gap mode is associated with q = 0.2027/d and E = 0.6682 × 10 −3 meV.

Odd N
The composition of the eigenvectors slightly changes for the odd case compared to the even N case Although both odd sized chains share the same spectrum, it is possible to find a linear combination of states which belongs to one chain only. The form of the extended states of the odd chains (Δ = ±t and E ± = 0) does not differ much from the one of the even chain and the entries of v α are x l where with k n d = nπ/(N + 1) (n = 1, . . . , N, n = (N + 1)/2). The exchange of a's and b's leads again to the coefficients for the chain β (see appendix C). The significant difference between even and odd N lies in the realization of the open boundary condition. Solving which leads to the momenta k n . An SSH-like chain with an odd number of sites hosts only a single zero energy mode, but α and β contribute each with one. We find on subchain α for Δ = ±t and on subchain β where l runs from 1 to (N + 1)/2. Regarding the statistics of the operators ψ † α , ψ † β associated to the states ψ α , ψ β , we proceed like for the even N case. The use of the SSH-like basis from equation (22) and the entries of the state ψ α yield now Again, the equations (67), (68) and (70) imply a perfect compensation of the A and B sublattice contributions, yielding (ψ † α ) 2 = 0 for E = 0. The zero energy mode, given by its entries in equation (71), leads to (ψ † α ) 2 = 1/2. Further, we find that both zero energy modes ψ α,β have their maximum at opposite ends of the Kitaev chain and decay into the chain. To better visualize this it is convenient to introduce the decay length and remembering that the atomic site index of x l is j = 2l − 1 equation (71) yields for |t| |Δ| For |t| |Δ| the x l coefficients are given by the same equation without the (−1) l−1 factor. We have moreover q = ±1/ξ, where q is the imaginary momentum yielding E = 0 in equation (41). Thus the localisation of these states is determined only by t and Δ. In the parameter setting of Δ = t we find: while both states exchange their position for

The particle-hole-operator
In the last section we have shown that some of the zero energy eigenstates of the BdG Hamiltonian of the finite Kitaev chain are Majorana zero modes (MZM) by exploiting the statistics of the corresponding operators ψ α,β . We further corroborate this statement now by recalling that an MZM is defined as an eigenstate of the Hamiltonian H and of the particle hole symmetry P. The latter acts on an eigenstate ψ α,β of energy E by turning it into an eigenstate of H of energy −E. Thus, the energy of such an exotic state has to be zero, since eigenstates associated to different energies are orthogonal.
The three symmetries, time reversal, chiral and the particlehole symmetry, discussed in section 2.3, can be constructed in real space too. Of particular interest is their representation in the SSH-like basis. The antiunitary particle-hole symmetry is The time reversal and the chiral symmetry depend on N. If N is even we find The expressions for odd N follow by removing the last line and last column in each diagonal block. The effect of P from equation (79) can be seen explicitly if one considers P ψ α . For N even and Δ = t the elements x l , y l of ψ α are given in equations (51) and (52). Here y l /x 1 is pure imaginary and x l /x 1 is real. Hence, ψ α is not an eigenstate of P since the prefactor to y l is finite, i.e., E ± = 0. We conclude that in a finite Kitaev chain with even number of sites Majorana zero modes emerge only at the Kitaev points for μ = 0, since the states in equations (61)-(64) are eigenstates of P as well. In the situation of odd N and μ = 0, the eigenstates given by their elements in equations (71), (72) are Majorana zero energy modes for an appropriate choice of x 1 . These states can be delocalised over the entire chain, depending on their decay length ξ, while the case of full localisation is only reached at the the Kitaev points, where the MZM turn into the states given by equations (75)-(78).

Spectrum
The last missing situation is to consider a finite chemical potential μ. For this purpose we use the so called chiral basiŝ where, however, h and h † do not commute except for t = 0 or Δ = 0. Thus, such matrices cannot be diagonalised simultaneously. Nevertheless the eigenvalues η j (η * j ) of h (h † ) are easily derived e.g. following [33]. We find independent of whether Δ t or t > Δ. 5 and we need only to focus on det(h). If a single eigenvalue η j of h is zero then det(h) vanishes. Thus, for a zero energy mode the chemical potential must satisfy Obviously, equation (87) cannot be satisfied for generic values of t 2 − Δ 2 , because all other quantities are real. The only possibility is t 2 Δ 2 . There is only one exception for odd N, because the value n = (N + 1)/2 leads to μ = 0 in equation (87) for all values of t and Δ, in agreement with our results of section 3. This result is exact and confirms findings from [21,22]; further it improves a similar but approximate condition on the chemical potential discussed by Zvyagin in [23]. An illustration of these discrete solutions μ n , which we dub 'Majorana lines', is shown in figure 8. All paths contain the Kitaev points at μ = 0 and t = ±Δ. Further, their density is larger close to the boundary of the topological phase, as a result of the slow changes of the cosine function around 0 and π.
For growing number of sites N, the density of solutions increases. In the limit N → ∞, θ n = nπ/N + 1 takes all values in [0, π] and the entire area between μ = ±2 t 2 − Δ 2 for t 2 Δ 2 is now occupied with zero energy modes.
Regarding the remaining part of the topological region, we are going to show in the next section that in that parameter space only boundary modes with finite energy exist. Because the energy of these modes is decreasing exponentially with the system size, in the thermodynamic limit their energy approaches zero and the full topological region supports zero energy modes.
where we defined v = (ξ 1 , ξ 2 , . . . , ξ N ) T . Notice that we are not really interested in the eigenvector v here; we simply use its entries as dummy variables to release a structure hidden in the product of h and h † . The elements of h, where n, m = 1, . . . , N, allow us to calculate the product hh † entry wise Thus, importantly, equation (88) reveals a recursion formula for the components of v. The entries ξ are a generalisation of the Fibonacci polynomials ζ j from equation (27), to which they reduce for μ = 0, and may be called Tetranacci polynomials [41,42]. Further, we find the open boundary conditions from equation (88) to be where we used equation (89) for simplifications. Appendix D contains the description of how to deal with those polynomials, the boundary conditions and further the connection of equation (89) to Kitaev's bulk spectrum λ = E ± (k) in equation (7). Essentially one has to use similar techniques as it was done for the Fibonacci polynomials, where now the power law ansatz ξ j ∝ r j leads to a characteristic equation for r of order four. Thus, we find in total four linearly independent fundamental solutions r ±1,±2 , which can be expressed in terms of two complex wavevectors denoted by κ 1,2 through the equality These wavevectors are not independent, but coupled via For μ = 0 we can recover from equation (92) our previous results, whereby one has only pure real (k) or pure imaginary (iq) wavevectors 5 . Further, equations (7) and (92) yield The linearity of the recursion formula equation (89) states that the superposition of all four fundamental solutions is the general form of ξ j . Since the boundary conditions translate into a homogeneous system of four coupled equations and a trivial solution for ξ j has to be avoided, we find that the determinant of the matrix describing these equations has to be zero. After some algebraic manipulations, this procedure leads finally to the full quantization rule of the Kitaev chain where we introduced the function F(κ 1 , κ 2 ) as Similar quantization conditions are known for an open X-Y spin chain in transverse field [26]. Notice that the quantization rule is symmetric with respect to κ 1,2 . Table 1 gives an overview of the quantization rules for different parameter settings (Δ, t, μ). The bulk eigenvalues of a finite Kitaev chain with four sites and μ = 0 are shown in figure 9. The previous relations open another route to finding the condition leading to modes with exact zero energy. A convenient form of equation (92) is and the dispersion relation can be transformed into Both combinations κ 1 ± κ 2 yield the same energy, due to equation (95). If one of the brackets in equation (96) vanishes for κ 1 + κ 2 (κ 1 − κ 2 ), the second one does so for κ 1 − κ 2 (κ 1 + κ 2 ) too. Hence, zero energy is achieved exactly if This puts restrictions on Δ, t, μ. Together with the quantization rule in equation (93), this ultimately leads to (87) and to the condition 2 t 2 − Δ 2 < |μ| < 2|t|, defining the region where exact zero modes can form.
Regarding the remaining part of the topological phase diagram, we find that in the limit N → ∞ the difference between even and odd N vanishes and the part of the μ = 0 axis between t = −Δ and t = Δ for even N leads to zero energy states too, in virtue of equation (42). Analogously, the area around the origin in figure 8, defined by 2 t 2 − Δ 2 < |μ| and μ < 2|t| with μ = 0 does not support zero energy modes for all finite N. Instead, this area contains solutions with exponentially small energies, see equation (96), which become zero exclusively in the limit N → ∞. The wavevectors obey κ 1 = iq 1 , κ 2 = π + iq 2 with real q 1,2 for Δ 2 > t 2 , and κ 1 = iq 1 , κ 2 = iq 2 otherwise, which follows from equation (93) after some manipulations; see also [26]. Thus, the entire non trivial phase hosts zero energy solutions for N → ∞.

Eigenstates
The calculation of an arbitrary wave function of the Kitaev Hamiltonian without any restriction on t, Δ, μ is performed at best in the chiral basis yielding the block off-diagonal structure in equation (82). A suitable starting point is to consider a vector w in the following form . . . , σ N ). Solving for an eigenstate with eigenvalue λ demands on v, u with h from equation (83). Thus, h h † v = λ 2 v, and we recover equation (88) and the entries of v obey equation (89) again. In appendix E we derived the closed formula for ξ j , namely where ξ −2 , . . . , ξ 1 are the initial values of the polynomial sequence dependent on the boundary conditions, and X i ( j) inherit the selective property That these functions X i ( j) exist and that they indeed satisfy equation (100) for arbitrary values of j is discussed in appendix E. The remaining task is to obtain the initial values ξ −2 , . . . , ξ 1 , which follow from the open boundary conditions equation (90). Further, one has one free degree of freedom, which we to choose to be the entry ξ 1 of v. In total our initial values are ξ 1 , ξ 0 = 0; ξ −1 = bξ 1 /a and ξ −2 follows from ξ N+1 = 0 and equation (100), Demanding further bξ N+2 − aξ N = 0 quantizes the momenta κ 1,2 and in turn λ = E ± (κ 1,2 ), according to equation (93) (the relation between X j and κ 1,2 is discussed in the appendix E). Notice that the form given by equations (100) and (102) (below) hold for all eigenstates of the Kitaev BdG Hamiltonian; the distinction between extended/decaying states and MZM is made by the values of k 1,2 , or equivalently κ 1,2 . The second part of the eigenstate w, i.e. u, follows in principle from equation (99). A simpler and faster way is to consider h † h u = λ 2 u. A comparison of h and h † reveals that they transform into each other by exchanging a and b and switching μ into −μ. Consequently, the structure of the entries of u follow essentially from the ones of v. Thus, we find that the σ i obey equation (89) too, yielding with the same functions X i (j). The boundary condition on u reads Proceeding as we did for v yields σ 0 = 0, σ −1 = aσ 1 /b and where σ 1 is fixed by the first line of equation (99) The last open boundary condition aσ N+2 − bσ N = 0 is satisfied for λ = E ± (κ 1,2 ) and thus already by the construction of v. We notice that we assumed λ = 0 in order to obtain σ 1 ; the case λ = 0 is discussed in section 4 below.
In the limit of μ = 0 on w it holds that for all values of l; thus only two initial values for ξ j and two for σ j are necessary to fix the sequences. Finally, one can prove easily that the functions X i (j) are always real and in consequence all ξ j (σ j ) are real (pure imaginary) if ξ 1 is chosen to be real. Thus, the corresponding operators ψ w (ψ † w ) can never square to 1/2.

MZM eigenvectors at finite μ
The technique outlined above (in particular, equation (104)) cannot be used directly for exact zero energy modes, because then λ = 0 in equations (98)   The zero energy values are twice degenerated, as one can see from equation (86), and the associated zero modes are connected by the chiral symmetry C. Thus, we get zero energy states by superposition ψ A,B := ψ ± C ψ /2 too. The chiral symmetry equation (80), contains an alternating pattern of ±1, such that ψ A ( ψ B ) includes only non zero entries on the Majorana sublattice A (B). Hence, ψ A ( ψ B ) contains only x l (X l ) and Y j (y j ) terms and the last component depends on whether N is odd or even. In the latter case we have The form of the odd N eigenvectors is quite similar, see equations (G15) and (G16).
The composition of ψ A is illustrated in figure 10, where its entries are shown to form a sawtooth like pattern, following the entries of ψ A on both SSH-like chains.
The full calculation is given in appendix G. We focus here on ψ A exclusively, because the ψ B components follow essentially from ψ A by exchanging a and b and replacing iμ by −iμ. The chemical potential has still to obey equation (87).
The components of the zero mode ψ A have to satisfy (l = 1, . . . , for even N, and the open boundary conditions are The situation for the entries of ψ A for odd N is similar Solving these recursive formulas leads in both cases to with θ n = nπ/(N + 1), n = 1, . . . , N and where x 1 is a free parameter and −a/b 0 due to t 2 Δ 2 .
Recalling that a = i (Δ − t), b = i (Δ + t), equations (109), (110) predict an oscillatory exponential decay of the coefficients x l , Y l . For example where the decay length is defined by ξ = d/ ln t−Δ t+Δ , for t > Δ > 0. Summarizing: the zero energy modes ψ A,B look like small or strongly suppressed standing waves with n − 1 nodes for n = 1, . . . , N max and N max = N/2 (N max = N − 1/2) for even (odd) N. The expressions for X l and y l are obtained in a similar way and X 1 can be freely chosen. The open boundary conditions for l = 0 are satisfied by construction of Y l (y l ), while the remaining ones follow due to the structure of θ n . The zero mode ψ A is shown in figure 11 for a various range of parameters. For not too large ratios t/Δ > 1, the zero mode ψ A is mostly localised at one end of the Kitaev chain and decays away from it in an oscillatory way. The eigenstate ψ B is concentrated on the opposite end. While the oscillation depends on the chemical potential μ n associated to the zero mode, according to equation (87), the decay length is only set by the parameters Δ and t. Thus, as the ratio of t/Δ is increased, the zero energy mode gets more and more delocalized.
The zero energy states ψ A,B are MZM's, since they are eigenstates of the particle hole operator P equation (79) for real or pure imaginary values of x 1 , X 1 . Further, the states ψ = ψ A + ψ B and C ψ = ψ A − ψ B are MZM's too. On the other hand a fermionic state is constructed with ψ ± = ψ A ± i ψ B , similar to what was found in section 4 for the μ = 0 case, or at the Kitaev points in equations (4a) and (4b) in section 2.1.
There are three limiting situations we would like to discuss: t → ±∞, N → ∞, and how the eigenstate changes if the sign of the chemical potential is reverted. For the first situation we notice that larger hopping amplitudes affect the decay length ξ. Because −a/b → 1 for t → ±∞, this implies also that ξ → ∞ in that limit. Hence oscillations are less suppressed for large values of t, as illustrated in figure 11. Already a ratio of t/Δ ≈ 100 is enough to avoid a visible decay for N ≈ 20. This  (109) and (110). The decay length of the MZM increases for larger ratios of t/Δ for a fixed value of θ n , until the state is delocalised over the entire system. Lowering the chemical potential, e.g. following the vertical orange line, but keeping t/Δ fixed, changes the shape of the MZM's. Large decay length and chemical potential leads to Majorana modes which have highest weight in the center of the chain. effect can be found as long as N is finite, for appropriate values of the ratio t/Δ.
What happens instead for larger system sizes? Regardless of how close −a/b is to 1, for a finite t, at some point the exponent j(l) in x l Y l , X l and y l leads to significantly large or small values. Thus, the state ψ A ( ψ B ) becomes more localised on the left (right) end for t > 0, and on the right (left) one for t < 0.
If we change the chemical potential to its negative value we find that y l , Y l only change their sign. For odd N and for θ n = π/2, i.e. μ = 0, one recovers the result in the equations (71) and (72).

Numerical results and impact of disorder
In this section we discuss the impact of disorder on the topological boundary states. To this extent we investigate numerically the lowest energy eigenvalues of the finite Kitaev chain.

The clean Kitaev chain
The features predicted analytically above are also clearly visible in the numerical calculations. The lowest positive energy eigenvalues E 0 of a finite Kitaev chain, with the Hamiltonian given by equation (1) and varying parameters, are shown in figure 12. The phase diagram in figure 12(a) is the numerical equivalent of that shown in figure 11, but for a smaller range of t and μ. Because of the necessarily discrete sampling of the parameter space, the zero energy lines are never met exactly, hence along the Majorana lines we see only a suppression of E 0 . Along the border of the topological regime, μ 2t, all the boundary states in a finite system have finite energy, as shown in figure 12(b). Figure 12(c) displays E 0 for fixed t/Δ = 17, as a function of μ and N. The number of near-zero energy solutions increases linearly with N, according to equation (87). It is worth noting that a spatial overlap between Majorana components of the end states in a short system does not need to lead to finite energy (cf. the right column of figure 11). The decay length ξ of the in-gap eigenstates, defined in equation (73), is determined by the ratio t/Δ and is  the same both for the near-zero energy states along the Majorana lines and for the finite energy states between them. It is the maximum energy of the boundary states that decreases as E 0,max ∝ exp(−Nd/ξ) as N is increased, as illustrated in figure 12(d), in agreement with equation (42) and [2,23,24]. The minimum energy of zero can be reached for any chain length, provided that the chemical potential is appropriately tuned.

Topological protection against Anderson disorder
One of the most sought after properties of topological states is their stability under perturbations which do not change the symmetry of the Hamiltonian. In order to see whether the non-Majorana topological modes enjoy greater or lesser topological protection than the true Majorana zero modes, we have calculated numerically the spectrum of a Kitaev chain with Anderson-type disorder as a function of μ for N = 20 and three different values of t/Δ. The disorder was modeled as an on-site energy term ε i whose value was taken randomly from the interval [−W, W]. The energies E 0 , E 1 of the two lowest lying states are plotted in figure 13. In all plots W = 4Δ, i.e. twice larger than the ±2Δ gap at μ = 0. Each curve is an average over 100 disorder configurations. Even with this high value of the disorder it is clear that the energy of the ingap states is rather robust under this perturbation. For t/Δ = 2, i.e. close to the Kitaev point where the boundary states are most localized, they are nearly immune to disorder-its influence is visible only at high μ and in the energy of the first extended state. For higher ratios of t/Δ, closer to the value of N + 1 (cf equation (40) and the discussion under equation (42)), the lowest energy states seem to be strongly perturbed and the Majorana zero modes entirely lost. This is however an artifact of the averaging-the energy E 0 plotted in figure 14 for several disorder strengths W shows that for any particular realization of the disorder the zero modes are always present, but their positions shift to different values of μ. The existence of these crossings is in fact protected against local perturbations and they correspond to switches of the fermionic parity [22].

Conclusion
Due to its apparent simplicity, the Kitaev chain is often used as the archetypal example for topological superconductivity in one dimension. Indeed, its bulk spectrum and the associated topological phase diagram are straightforward to calculate, and the presence of Majorana zero modes (MZM) at special points of the topological phase diagram, known as Kitaev points (μ = 0, t/Δ = ±1 in the notation of this paper), is easy to demonstrate. However, matters become soon complicated when generic values of the three parameters μ, Δ and t are considered.
In this work we have provided exact analytical results for the eigenvalues and eigenvectors of a finite Kitaev chain valid for any system size. Such knowledge has enabled us to gain novel insight into the properties of these eigenstates, e.g. their precise composition in terms of Majorana operators and their spatial profile.
Our analysis confirms the prediction of Kao [21], whereby for finite chemical potential (μ = 0) zero energy states only exists for discrete sets of μ(Δ, t) which we dubbed 'Majorana lines'. We calculated the associated eigenvectors and demonstrated that such states are indeed MZM. Importantly, such MZM come in pairs, and because they are made up of Majorana operators of different types, they are orthogonal. In other words the energy of these modes is exactly zero, even when the two MZM are delocalized along the whole chain (which depends on the state's decay length ξ).
Beside of the Majorana lines, but still inside the topological region, finite energy boundary states exist. We studied the behavior of the energy E 0 of the lowest state numerically as a function of t/Δ, μ/Δ and of the system size L = Nd. We found that with good accuracy the maximum energy E 0,max ∝ exp(−L/ξ). This energy, and hence the energy of all the boundary states, tends to zero in the thermodynamic limit N → ∞. For fixed N the ratio t/Δ can be varied until the decay length ξ becomes of the order of the system size L and hence the associated E 0,max is not exponentially close to zero energy.
All the boundary states in the topological region, whether of exact zero energy or not, are of topological nature, as predicted by the bulk-edge correspondence. This fact is important in the context of topological quantum computation. In fact, whether a state has exact zero energy or not is not relevant for computation purposes, as long as this state is topologically protected. For the prototypical model of the finite Kitaev chain considered in this paper this condition is best satisfied in the vicinity of the Kitaev points t/Δ = ±1, μ = 0 independent of the system's size. When a random disorder of Anderson type is introduced in the Kitaev chain, our numerical results show that the upper bound E 0,max is remarkably robust, confirming the topological protection of all topological states.
The fact that even in short chains, with large spatial overlap of the Majorana components, the energy of the subgap state can still be strictly zero could have interesting implications for the experiments. For example, in STM-spectroscopy experiments with topological chains like in Nadj-Perge et al in [12] or Kim et al in [16], this would yield a zero-bias Majorana peak delocalized along the whole chain.
Although our treatment using Tetranacci polynomials for μ = 0 is general (sections 5 and 6), we have dedicated special attention to two parameter choices in which the Tetranacci polynomials reduce to the generalized Fibonacci polynomials. The first case is that of zero chemical potential discussed in sections 3 and 4 of the paper, where the Kitaev chain turns out to be composed of two independent SSH-like chains. This knowledge allows one a better understanding of the difference between an even and an odd number N of sites of the chain. This ranges from different quantization conditions for the allowed momenta of the bulk states, to the presence of MZM. While MZM are always present for odd chains, they only occur at the Kitaev points for even chains. When μ is allowed to be finite, the Kitaev points develop into Majorana lines hosting MZM for both even and odd chains. In the thermodynamic limit the distinction between even and odd number of sites disappears. The fact that E = 0 at the Majorana lines decouples the two Majorana sublattices, and again allows us to use the simpler Fibonacci polynomials.
Deutsche Forschungsgemeinschaft via SFB 1277 Project B04. We acknowledge useful discussions with A Donarini, C de Morais Smith and M Wimmer.

Appendix A. A note on Fibonacci and Tetranacci polynomials
An object of mathematical studies are the Fibonacci numbers F n (n ∈ N 0 ) defined by which frequently appear in nature. A more advanced sequence is the one of the Fibonacci polynomials [36] F n (x), where with an arbitrary complex number x which gives different weight to both terms. The polynomial character becomes obvious after a look at the first few terms The so called generalized Fibonacci polynomials [37] are defined by where x, y are two complex numbers. The second weight changes the first elements of the sequence in equation (A2) to 0, 1, x, x 2 + y, x 3 + 2x y.
There is a general mapping between the sequences F n (x) in equation (A2) and F n (x, y) in equation (A4), namely where F n (x/ √ y) obeys equation (A2) with x/ √ y instead of x. A last generalization is to consider arbitrary initial values F i (x, y) = f i , i = 0, 1 and keeping [38] This changes the first terms into , The Fibonacci polynomials ζ 2n−1 , ζ 2n , 2n−1 , 2n we consider in the spectral analysis are of the last kind with x = λ 2 + a 2 + b 2 , y = −a 2 b 2 (a = i(Δ − t), b = i(Δ + t)) and with different initial values for odd and even index as well as different ones for ζ and . The first terms for ζ 2n−1 are while one finds for ζ 2n ζ 0 = 1, The expressions for 2n ( 2n−1 ) follow from the ones of ζ 2n (ζ 2n−1 ) by exchanging a and b. A closed form for Fibonacci numbers/polynomials is called a Binet form, see for example [38]. In the case of ζ n this form is given in equations (B20) and (B21).
In order to obtain the general quantization condition for the wavevectors of the Kitaev chain we face further generalizations of Fibonacci polynomials, so called Tetranacci polynomials τ n , defined by with four complex variables x 0 , . . . , x 3 and four starting values τ 0 , . . . , τ 3 . These polynomials are a generalisation of Tetranacci numbers [41,42] and their name originates from the four terms on the r.h.s. of equation (A8). The form of Tetranacci polynomials ξ j we deal with in this work, is provided by equation (D4).

B.1. Characteristic polynomial in closed form
The full analytic calculation of the spectrum is at best performed in the basis of Majorana operators γ A(B) j , ordered according to the chain index Then the BdG Hamiltonian becomes block tridiagonal where A and B are 2 × 2 matrices Since we are interested in the spectrum, we have essentially only to calculate (and factorise) the characteristic polynomial P λ H KC = det (λ − H KC ) which reads simply (B4) at zero μ. In the following we will consider λ to be just a 'parameter', which is not necessarily real in the beginning. Further, we shall impose (only in the beginning) the restriction λ = 0. However, our results will even hold without them. The validity of our argument follows from the fact that the determinant P λ is a smooth function in the entire parameter space P := C 3 , which contains a, b and λ.
The technique we want to use to evaluate P λ is essentially given by the recursion formula of the 2 × 2 matrices Λ j [43,44]: where j = 1, . . . , N and P λ = N j=1 det Λ j . The matrices B and B † are pure off-diagonal matrices and since λ 2 is diagonal, one can prove that Λ j has the general diagonal form of Λ j := x j 0 0 y j (for all j). The application of equation (B6) leads to a recursion formula for both sequences of entries and the initial values are x 1 = y 1 = λ. We find x j and y j to be fractions in general, and define ζ j , j , β j and δ j by to take this into account. The initial values can be set as and after a little bit of algebra we find their growing rules to be where j starts from 1. The definitions ζ 0 := δ 1 = 1 and 0 := β 1 = 1, enable us to get rid of the δ j and β j terms inside equations (B9) and (B10). Hence which leads to the relations We already extended the sequences of ζ j and j artificially backwards and we continue to do so, using the equations (B13) and (B14), starting from j = −1 with ζ −1 = −1 = 0. Please note there are no corresponding x 0 , y 0 or even x −1 , y −1 expressions, since they would involve division by 0.
The last duty of β j and δ j is to simplify the determinant P λ by using the equations (B8), (B11) and (B12) which reduces the problem to finding only ζ N and N . Please note that the determinant is in fact independent of the choice of the initial values for ζ 1 , 1 , β 1 and δ 1 in the equations (B7) and (B8). Further, equations (B13), (B14) and (B17) together show the predicted smoothness of P λ in P and all earlier restrictions are not important anymore. Finally we consider λ to be real again.
Even though it seems that we are left with the calculation of two polynomials, we need in fact only one, because both are linked via the exchange of a and b. Note that λ is considered here as a number and thus does not depend on a and b. Further, the dispersion relation is invariant under this exchange.
The connection of ζ j and j for all j −1 is and can be proven via induction using equations (B13) and (B14). Decoupling ζ j and j yields where one identifies them as (generalized) Fibonacci polynomials [36,37]. The qualitative difference between even and odd number of sites is a consequence of equation (B18) and the initial values for ζ j . The next step is to obtain the closed form expression of ζ j ( j ), the so called Binet form. We focus exclusively on ζ j .
The Binet form can be obtained by using a power law ansatz u j ∝ r j , leading to two fundamental solutions Please note that this square root is always well defined, which can be seen in the simplest way by setting λ to zero. Consequently, the difference between r 1 and r 2 is never zero. A general solution of u j (v j ) can be achieved with a superposition of r 1,2 with some coefficients c 1,2 , due to the linearity of their recursion formula. Both constants c 1 and c 2 are fixed by the initial values for ζ j , for example u 0 = ζ −1 = 0 and u 1 = ζ 1 = λ and similar for v j . After some simplifications, we finally arrive at in agreement with [36][37][38]. The validity of the solutions is guaranteed by a proof via induction, where one needs mostly the properties of r 1,2 to be the fundamental solutions. The exchange of a and b leads to the expressions where we used that r 1,2 is symmetric in a and b. At this stage we have the characteristic polynomial in closed form for all Δ, t and more importantly for all sizes N at zero μ. We can already anticipate the twice degenerated eigenvalues of the odd sized Kitaev chain, because from the closed forms of j and ζ j it follows immediately Notice that equation (B24) is important to derive the characteristic polynomial via the SSH description of the Kitaev BdG Hamiltonian at μ = 0 and to show the equivalence to the approach used here. It is recommended to use the determinant formula in [45] together with equations (B9) and (B10) for the proof.
The main steps of the factorisation are mentioned in the next section.

B.2. Factorisation of generalized Fibonacci polynomials
The trick to factorise our Fibonacci polynomials [36,37] bases on the special form of r 1,2 . The ansatz is to look for the eigenvalues λ in the following form which is actually the definition of θ. The hermiticity of the Hamiltonian enforces real eigenvalues and consequently θ can be chosen either real, describing extended solutions, or pure imaginary, which is connected to decaying states. The ansatz leads to an exponential form of the fundamental solutions and we consider θ ∈ R first. Thus, we find the eigenvalues for odd N One obvious solution is λ = 0. The introduction of 2kd = θ, where d is the lattice constant of the Kitaev chain, leads to: and solutions inside the first Brillouin-zone are given by where n runs from 1, . . . , N without (N + 1)/2. Please note that equation (B26) cannot be satisfied for N = 1.
The even N case requires more manipulations. We first rearrange equation (B21) as The expressions λ 2 + b 2 − r 1,2 are simplified to In the end ζ 2j becomes Note that the competition of Δ and t is hidden inside the square root affecting both the quantization condition and the dispersion relation E ± (k) = λ(θ), which follow from equation (B25). However, both situations lead to the same result, because the momenta and the spectrum are shifted by π/2 (with respect to kd). From ζ N it follows:  (34) and (35). The case of decaying states is similar, but not just done by replacing k by iq. The following case is only valid for even N, since we have already all 2N eigenvalues of the odd N case.
Our ansatz is modified to by an additional minus sign, which is important to find the decaying state solutions. After some manipulations ζ N = 0 yields the quantization conditions where qd = θ/2. The conditions for qd = 0 as solution, corresponding to infinite decay length ξ, turn out to be ±t/Δ = N + 1 (if |t| |Δ|) or ±Δ/t = N + 1 (else) and follow by applying the limit qd → 0 on equations (B30) and (B31). A last simplification can be done for qd = 0 The criterion to find a wave vector is that (−m) 1, but not larger than N + 1, which leads then to exactly two solutions ±q and otherwise to none. The corresponding eigenvalues can be obtained from E ± (qd) = ± 4t 2 cosh 2 (qd) − 4Δ 2 sinh 2 (qd), |Δ| |t|, E ± (qd) = ± 4Δ 2 cosh 2 (qd) − 4t 2 sinh 2 (qd), |t| |Δ|, which can be zero. The results for N can be obtained by replacing t with −t everywhere.
Due to the linearity of the recursion formula, the most generic ansatz for y l is y l = c 1 f l 1 + c 2 f l 2 , where the constants c 1,2 follow from the initial values of y 1,2 . The calculation of both constants leads to where T l is simply [40] T l : Analogously we find A short comment on the initial values y 1,2 . A hermitian matrix is always diagonalisable, regardless of degeneracies in its spectrum and an eigenvector is well defined only up to the prefactor. Consequently we have the freedom to choose one component of v α . This choice will in turn define all remaining initial values.
Consider for example x 1 to be a fixed value of our choice. We find y 1 , x 2 and y 2 to be The y 2 can be rewritten as y 2 = y 1 [ f 1 + f 2 ] which leads to a simpler form of all y l 's [40] After a bit of algebra, one finds x l to be So far we found the general solutions of the recursion formulas equations (C1)-(C4). The comparison of equations (C2) and (C3) leads to because the recursion formulas themselves do not care about any index limitation. The last equation means only that the wave function of a finite system has to vanish outside, at the boundary, yielding the quantization rule.
The extended states can be obtained with where equations (C6) and (C7) relate kd and λ, and T l is recast as The last equation for T l yields via equations (C10) and (C11) the quantization condition. Thus, the momenta k obey where each solution defines two states with the energy E ± from equation (35).
The decaying states depend strongly on the interplay of Δ and t. The ansatz is where s is defined as Finally the coefficient T l becomes The proper q, if existent, leads to two states and satisfies where m is In total we have already all N non normalized states with respect to the chain α and this approach holds as long as |Δ| = |t|. The remaining cases start again from the equations (C1)-(C4). Case 2. Eigenvectors at the Kitaev point. We consider now Δ = −t, or b = 0, and we have to solve where l runs from 1 to (N − 2)/2. A zero energy mode is obviously not existing on the α subchain, because λ = 0 would lead to v α = 0 which is not an eigenvector by definition. These zero modes belong to the subchain β for Δ = −t. The only possible eigenvalues for the extended modes of the α chain are λ = ±2t [1,2], see equation (35). Recalling a = −2it, leads to N/2 independent solutions of dimerised pairs (x l , y l ) with y l = ∓ix l and the signs are with respect to the eigenvalues.
The last cases belong to Δ = t (a = 0), where we search for the solution of where l runs from 1 to (N − 2)/2. The first (second) line clearly states that either λ is zero and/or x 1 (y N/2 ). The zero λ means on the one hand that most entries vanish x 2 = x 3 = · · · = x N/2 = 0 and y 1 = y 2 = · · · = y (N−2)/2 = 0, since b = 2it = 0 to avoid a trivial Hamiltonian. On the other hand we have two independent solutions, first and second describing the isolated MZM's at opposite ends of the chain. In the case of λ = ±2t, we have N − 2 independent solutions in form of pairs (y l , x l+1 ) with y l = ±ix l+1 .
The non trivial solutions for v β follow by replacing x l → X l , y l → Y l and t → −t everywhere.

C.2. N odd
The eigenvectors have similar shape v α = x 1 , y 1 , x 2 , y 2 , . . . , , but the last entry is different compared to the even N case. Although both subchains have the same spectrum, it is possible to consider a superposition of eigenstates of the full Hamiltonian which belongs to only one chain, for example α. We consider v β to be zero.
The eigenvector system for v α reads with l = 1, . . . , N−3 2 and i = 1, . . . , N−1 2 . If we consider a, b and λ all to be different from zero, we find again that the entries of v α are Fibonacci polynomials obeying the same recursion formula as in the even N case and lead to the same solution where l = 1, . . . , N+1 2 and T l,(i) is as before. The ansatz f 1 = e 2ikd , f 2 = e −2ikd for the extended states influences T l (T i analogously) and leads via y N+1 to the equidistant quantization k ≡ k n = n π N+1 with n = 1, . . . , (N − 1)/2, due to the number of eigenvectors of the single SSH-like chain. Both chains share the same spectrum for odd N and thus we have in total n = 1, . . . , N, n = (N + 1)/2.
We report here shortly on all other parameter situations.
(a) If we consider a and b to be different from zero, but λ = 0, we find only one state and l runs from 1 to (N − 1)/2. The results for the α chain follow again by replacing x l → X l , y l → Y l and t → −t.

Appendix D. Spectrum for finite μ
The BdG Hamiltonian, expressed in the chiral basisΨ As mentioned in section 5, we look for a solution of with v = (ξ 1 , ξ 2 . . . , ξ N ) T to find the general quantization rule. The entries of the matrix hh † are and equation (D3) becomes the Tetranacci sequence where j = 1, . . . , N − 5. The missing four boundary terms are We extend the Tetranacci sequence from j = −∞ to j = ∞, i.e. the index limitations in equation (D4) can be ignored, while v still contains only ξ 1 , . . . , ξ N . Consequently, we can simplify the boundary conditions by using the recursion formula and further any restriction like N > 3 does not exist. We find The procedure we followed in the context of Fibonacci polynomials was to obtain a closed form with the ansatz ξ j = r j , r = 0. So we do here on starting from equation (D4). Thus, the characteristic equation for r reads and we have to find all four zeros to determine ξ j in the end. We introduce two new variables to simplify the expressions in the following. The characteristic equation becomes Dividing by r 2 (r = 0) and defining S := r + r −1 leads to where we can read out the solutions S 1,2 The definition of S amounts to an equation for r Thus one can insert the solutions S 1,2 and solve for r. We find yielding directly r j r −j = 1. Here, we choose the ansatz which is actually the definition of κ 1,2 . Since the coefficients S 1,2 contain λ through the variable ζ, this is in the end an ansatz for λ. The expression for λ follows easily from equation (D10) by inserting equation (D13). Using the definiton of η and resolving for ζ first and in a second step for λ we finally arrive at Kitaev's bulk formula λ κ 1,2 = ± μ + 2t cos κ 1,2 2 + 4Δ 2 sin 2 κ 1,2 .
Let us return to ξ j . Since the recursion formula in equation (D4) is linear, a superposition of all four solutions r ±1 , r ±2 is still a solution with some coefficients c 1,2,3,4 ∈ C. From equation (D12) it follows and thus ξ n = c 1 e iκ 1 n + c 2 e −iκ 1 n + c 3 e iκ 2 n + c 4 e −iκ 2 n .
Further, equation (D15) implies that we consider a combination of states of the same energy. The usually following step would be to fix these constants, requiring four initial values. We can use e.g. ξ 1 as free parameter. Further setting ξ 0 = ξ N+1 = 0, ξ −1 = (b/a)ξ 1 as the boundary conditions yield a sufficient number of constraints.
The remaining condition aξ N = bξ N+2 yields the quantization rule then. However, if one is not interested in the state v or in the general eigenstates of the Kitaev chain, but only in the quantization rule, one can use a much simpler approach. Using our ansatz for ξ j from equation (D18) and being aware of the fact that the boundary conditions yield a homogeneous system, we find Demanding det (B) = 0 avoids a trivial solution and leads to the quantization rule in equations (93) and (94).

Appendix E. The closed formula of Tetranacci polynomials
The goal here is to obtain the general solutions of a polynomial sequence ξ j , j ∈ Z, which obeys with arbitrary initial values. We consider here ξ −2 , . . . , ξ 1 , other choices are possible too, to be the initial values. We want to determine a closed form expression for all ξ j 's. Similar to equation (D18), the general solution is given by a superposition of the four fundamental solutions r ±i from equation with some constants c 1 , . . . , c 4 which follow from ξ −2 , . . . , Solving equation (E3) and factorising ξ j into contributions of ξ −2 , . . . , ξ 1 yields The functions X i (j) obey by construction (via equation (E3)) for all values of λ, μ, a, b and further obey equation (E1).
Despite the short form of ξ j in equation (E4), the formulas of X i (j) tend to be lengthy, such that we first introduce a short hand notation for their main pieces. We define where the r.h.s of both equalities arise due to r i r −i = 1 for i = 1, 2. With S 1,2 from equation (D11) we find the X i ( j) to be whereσ is meant as 'not σ', e.g. if σ = 1 then we haveσ = 2 and vice versa. The presence of S 1,2 in the form of X i (j) arises due to the definition of F 1,2 , since they are Fibonacci polynomials with inital values F 1,2 (0) = 0, F 1,2 (1) = 1 and obey The proof is done by induction over j and using the relation between r ±i and S i according to the equations (D11) and (D12). The formulas for X i ( j) and F i are exact and hold for all values of μ, a, b (t, Δ) and for all values of λ, regardless of whether equation (E1) describes an eigenvector problem or not. Notice that μ = 0 is a special situation, since for μ = 0. Thus, we find for all values of l.
The closed formula of ξ j can be used in multiple ways. In the context of eigenvectors the exponential form of the fundamental solutions r ±i according to equation (D17) is the direct connection to the momenta κ 1,2 , and their values follow from the quantisation rule in equation (93). The corresponding value of λ = E ± (κ 1,2 ) follows then from equation (7). The form of F 1,2 transforms into a ratio of sin(κ 1,2 j)/sin(κ 1,2 ). However, once the energy E ± (κ 1,2 ) is known, the explicit use of κ 1,2 is not important, since r ±i follow also directly from equation (D12).

Appendix F. The zeros of the determinant
We notice the Fibonacci character [36][37][38] of the sequence in equation (F2) and continue with the calculation of the Binet form. The ansatz h j ∝ R j (R ∈ C\ {0}) leads to R 2 + iμ R − ab = 0, and the solutions R 1,2 are Our ansatz holds for all parameter choices of μ, Δ and t and R 1,2 obey R 1 + R 2 = −iμ, ( F 5 ) The general form of h j is given by a superposition of R 1 and R 2 h j = n 1 R j 1 + n 2 R j 2 , ( F 7 ) and n 1,2 are fixed by the initial values. The calculation can be simplified by extending the sequence h j backwards with equation (F2), because h −1 = 0. The use of h −1 and h 0 leads to yielding the closed form of h j We find the determinant of the Kitaev chain to be for all values of μ, t, Δ ∈ R. The determinant does not vanish in general, due to equation (F4), but only for a specific combination of the parameters μ, t, Δ.
In the following we consider t and Δ to be fixed values of our choice and we search for the values of μ such that the determinant vanishes. The Fibonacci character of h N enables us to factorize the determinant [36,37] and leads automatically to the zeros. The factorization follows from equation (F4) and the starting point is the square root: 4 ab − μ 2 = 4 (t 2 − Δ 2 ) − μ 2 .
We found all zeros in case (a) and we continue with (b).

F.3. Case (c)
We consider here Δ 2 t 2 and 4(t 2 − Δ 2 ) − μ 2 0. We start by manipulating the square root in R 1,2 4 ab − μ 2 = 4 (t 2 − Δ 2 ) − μ 2 = i μ 2 + 4 (Δ 2 − t 2 ). (F16) Our ansatz is μ = 2 Δ 2 − t 2 v(θ) with a real valued function v(θ), without further restrictions, because μ 2 = 4 (Δ 2 − t 2 ) v 2 (θ) −4 (Δ 2 − t 2 ), in view of μ 2 −4(Δ 2 − t 2 ). The square root in R 1,2 becomes in general i μ 2 + 4 (Δ 2 − t 2 ) = i (Δ 2 − t 2 ) v 2 (θ) + 1, and one sees immediately that v(θ) = sinh(θ), θ ∈ R is an appropriate choice. We find for R 1,2 the form where the negative sign in front of the exponential forces us to distinguish between even and odd N.  Consequently we need first of all |R 1 | = |R 2 |. The second part is to find the proper phase factors and all of them lie on a circle with radius |R 1 | in the complex plane. We have found non trivial solutions only for scenario (a). In total, we found all conditions det (H KC ) = 0. The general case is when the chemical potential is = 0. A last check for the odd N case is done by choosing n = (N + 1)/2, i.e. θ n = π/2, which leads back to the old μ = 0 limit. Applying θ n → π/2 on x l leads to after some steps, while all Y l ∝ μ are zero. Similar we find X l from x l upon changing a with b, while y l = 0 for all l. Hence, we recover our result for the α (β) chain, see equations (71) and (72). The remaining questions is whether these zero energy modes are Majorana zero modes or not. The use of the particle hole operator in the SSH-like basis from equation (79), i.e. complex conjugation, reveals that the expressions x l /x 1 , Y l /x 1 , X l /X 1 and y l /X 1 are always real quantities, for both even and odd N. Thus ψ A ( ψ B ) is an MZM if x 1 (X 1 ) is either real or pure imaginary.
The MZM mode ψ A ( ψ B ) has non-zero weight on γ A j (γ B j ). Superpositions of both vectors can be MZM too if the coefficients are chosen properly. For example ψ = ψ A + ψ B has no zero entry. Hence, it is a mixed type MZM (for the correct choice of x 1 and X 1 ).