\section{Introduction}
The quest for topological quantum computation has drawn a lot of attention to Majorana zero energy modes (MZM), quasi-particles obeying non-Abelian statistics hosted by topological superconductors \cite{Aguado}.
The archetypal model of a topological superconductor in one dimension was proposed by Kitaev \cite{kitaev:physusp2001}. It consists of a chain of spinless electrons with nearest neighbor 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 \cite{kitaev:physusp2001}.
This model has thus become very popular due to its apparent simplicity and it is often used to introduce topological superconductivity in one dimension \cite{Alicea-2010, Aguado}. Also more sophisticated realizations of  effective p-wave superconductors, based on semiconducting nanowire-superconductor nanostructures \cite{Lutchyn-2010, Oreg-2010, Mourik-2012, Jelena-2012, Deng-2016, Jelena-2017, Zhang-2018,Prada2019},
 ferromagnetic chains on superconductors\cite{Nadj-Perge-2013, Nadj-Perge-2014,Jelena-2013, Zvyagin-2013, Kim-2018}
 or s-wave proximitized carbon nanotubes \cite{Sau-2013, Marganska-2018, Milz-2018}, 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 Ref.~[\onlinecite{Nadj-Perge-2014}] it is expected that the chain length is shorter than the superconducting coherence length \cite{Lee-2014}.
Spectral properties of a finite-length Kitaev chain have been addressed in more recent papers  \cite{Kao-2014,Hegde-2015,Zvyagin-2015,Zeng-2019}, 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 \cite{Kao-2014}, a finite-length Kitaev chain also supports modes with exact zero-energy.
 However, they are only found for discrete values of the chemical potential. As shown by Hegde et al., in Ref.~[\onlinecite{Hegde-2015}], 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 \cite{Zvyagin-2015} 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 \cite{Lieb-1961, Loginov-1997}. Kawabata et al. in Ref.~[\onlinecite{Kawabata-2017}] 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  {\it{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 Refs. [\onlinecite{Hegde-2015}], [\onlinecite{Zvyagin-2015}] and respectively [\onlinecite{Zeng-2019}]. They also complete the analysis of the eigenstates of an open Kitaev chain performed by Kawabata et al.\cite{Kawabata-2017} 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 \cite{Aguado}. 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 $\Delta$, 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 (self-conjugate) 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\cite{kitaev:physusp2001} or the winding number\cite{Review-Chiu-2016} topological invariant. For a finite-length chain MZM only exist for a set of discrete values of the chain parameters, see Eq. (\ref{equation: eigenvalue zero condition kouachi}) below, in line with Refs. [\onlinecite{Kao-2014}] and [\onlinecite{Kawabata-2017}]. 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\cite{Kempkes-2019}. \\

The paper is organized as follows. Section \ref{section: The Kitaev chain and its bulk properties} shortly reviews the model and its bulk properties. Section \ref{section: spectral analysis} 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 (SSH-like) chains. \cite{wakatsuki:prb2014,Kawabata-2017} 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 \ref{section: Eigenvectors and symmetries}.  In Section \ref{section: influence of non zero chemical potential}, \ref{section: MZM paths in the topological phase diagram for finite N} and \ref{section:  numerical results and impact of disorder} we turn to the general case of finite chemical potential which couples the two SSH-like chains. While Sec. \ref{section: influence of non zero chemical potential} deals with the energy eigenvalues and eigenvectors of the finite chain, Sec. \ref{section: MZM paths in the topological phase diagram for finite N} provides exact analytical results for the MZM. In Sec. \ref{section:  numerical results and impact of disorder} the influence of disorder on the energy of the lowest lying states is investigated numerically. In Sec. \ref{section: conclusion} conclusions are drawn. Finally, appendices \ref{appendix: comments on Fibonacci and Tetrannaci polynomials}-\ref{appendix: zero energy states} contain details of the factorisation of the characteristic polynomial in real space and the calculation of the associated eigenstates.

\section{The Kitaev chain and its bulk properties}
\label{section: The Kitaev chain and its bulk properties}

\subsection{Model}
\label{section 2 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 $\mu$, the hopping amplitude $t$, and the p-wave superconducting pairing constant $\Delta$. The Kitaev Hamiltonian, written in a set of standard fermionic operators $d_j,\,d_j^\dagger$, is \cite{kitaev:physusp2001, Aguado}
\begin{align}\label{equation: Kit. Hamiltonian, fermionic operators, realspace}
	\hat{H}_{\mathrm{KC}}=-\mu\sum\limits_{j=1}^N\,d_j^\dagger d_j+\sum\limits_{j=1}^{N-1}\,\left(\Delta\, d_j^\dagger d_{j+1}^\dagger\,-t\,d^				\dagger_{j+1} d_{j}+h.c.\right),
\end{align}
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 $\Delta$ and $t$ to be real parameters from now on.

The Hamiltonian in Eq. (\ref{equation: Kit. Hamiltonian, fermionic operators, realspace}) 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\cite{kitaev:physusp2001}. This can be seen by expressing the Kitaev Hamiltonian in terms of so called Majorana operators $\gamma^{A,B},$
\begin{align}\label{equation: Definition Majorana operators}
	\left(\begin{matrix}
	d_j\\
	d_j^\dagger\end{matrix}\right)
	\,\mathrel={\mathop:}\,
	\frac{1}{\sqrt{2}}\,
	\left(\begin{matrix}
	1& i\\
	1 & -i
	\end{matrix}\right)\,
	\left(
	\begin{matrix}
	\gamma_j^A\\
	\gamma_j^B\end{matrix}
	\right), \quad \left(\gamma^{A,B}\right)^\dagger = \gamma^{A,B},
\end{align}
yielding the form
\begin{align}\label{equation: Kitaev chain in Majorana operators}
	\hat{H}_\mathrm{KC}\,&=\,-i\,\mu\,\sum\limits_{j=1}^N\,\gamma_j^A\gamma_j^B\,+\,i\left(\Delta+t\right)\sum\limits_{j=1}^{N-1}				\gamma_j^B\gamma_{j+1}^A \notag\\
	&\qquad +\,i\left(\Delta-t\right)\sum\limits_{j=1}^{N-1} \gamma_{j}^A \gamma_{j+1}^B.
\end{align}
Notice that, in virtue of Eq. \eqref{equation: Definition Majorana operators} it holds $\{\gamma_j^A,\,(\gamma_j^A)^\dagger\}\,=\,2\,(\gamma_j^A)^2\,=\,1$, and similarly for $\gamma_j^B$. For the particular parameter settings $\Delta = \pm t$ and $\mu=0$, which we call the Kitaev points, Eq. (\ref{equation: Kitaev chain in Majorana operators}) leads to a "missing" fermionic quasiparticle $q_\pm$:
\begin{subequations}
	\begin{align}
		\label{equation: Missing electron a}
		q_+\,&=\,\frac{1}{\sqrt{2}}\left(\gamma_1^A\,+\,i\,\gamma_N^B\right)\qquad \left[\Delta=t\right],\\
		\label{equation: Missing electron b}
		q_-\,&=\,\frac{1}{\sqrt{2}}\left(\gamma_1^B\,+\,i\,\gamma_N^A\right)\qquad \left[\Delta=-t\right].
	\end{align}
\end{subequations}
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 ($\mu=0$, $\Delta=\pm 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 $\mathcal{P}$. Equivalently, if $\gamma_M^\dagger$ is an operator associated to such a MZM it satisfies $\gamma_M^\dagger  =\gamma_M$ and $\gamma_M^2=1/2$.

The topological phase diagram is shortly reviewed in Sec. \ref{subsection: Topological phase diagram}.
\subsection{Bulk spectrum}
The Hamiltonian from Eq. (\ref{equation: Kit. Hamiltonian, fermionic operators, realspace}) in the limit of $N\rightarrow \infty$ reads in $k$ space
\begin{align}\label{equation: Kitaev Hamilt. in electronic operators}
	\hat{H}_\mathrm{KC}\,=\,\frac{1}{2}\sum_k\,\hat{\psi}^\dagger_k\,\mathcal{H}(k)\hat{\psi}_k,~~
	\hat{\psi}_k\,=\,\left(d_k,\,d_{-k}^\dagger\right)^\mathrm{T},
\end{align}
where we introduced the operators \makebox{$d_k=\frac{1}{\sqrt{N}}\sum\limits_j e^{-i\,j\,kd}\,d_j$} and $k$ lies inside the first Brillouin zone, i.e. \makebox{$k\in\left[-\frac{\pi}{d},\frac{\pi}{d}\right]$} and $d$ is the lattice constant. The $2\times 2$ Bogoliubov- de Gennes (BdG) matrix
\begin{align}\label{eq:hamiltonian-Kitaev-k}
	\mathcal{H}(k)\,=\,\left[\begin{matrix}
	-\mu-2t\cos(kd) & -2i\Delta\sin(kd)\\
	2i\Delta\sin(kd) & \mu+2t\cos(kd)
	\end{matrix}\right]
\end{align}
is easily diagonalized thus yielding the excitation spectrum
\begin{align}\label{infinite chain: excitation spectrum}
	E_\pm(k)\,=\,\pm\sqrt{4\Delta^2\,\sin^2(kd)\,+\,\left[\mu\,+\,2t\,\cos(kd)\right]^2}.
\end{align}
Note that for $\mu=0$ Eq. (\ref{infinite chain: excitation spectrum}) predicts a gapped spectrum whose width is either $4\,\Delta$ ($\vert \Delta \vert < \vert t \vert$) or $4\,t$ ($\vert t \vert < \vert \Delta \vert$).
\subsection{Topological phase diagram}
\label{subsection: Topological phase diagram}
\begin{figure}[h]

 \includegraphics[width=\columnwidth]{figure1.pdf}
	\caption{Topological phase diagram of the Kitaev chain for $\Delta>0$, constructed with the winding number invariant Eq. (\ref{equation: winding number invariant}). Distinct topological phases are separated by the phase boundary at $\mu = \pm 2t$, visualised by the red lines. The four insets illustrate the evolution of $w(k)$ along the Brillouin zone at the four marked positions in the phase space, ($t/\Delta =\pm 2$, $\mu/\Delta=0.2$) in the topological and ($t/\Delta =\pm 2$, $\mu/\Delta=4.2$) in the trivial phase. In the phase diagram for $\Delta<0$ the $\nu=+1$ and $\nu=-1$ regions are swapped.}
	\label{fig:topo-phase-diagram}
\end{figure}
The BdG Hamiltonian \eqref{eq:hamiltonian-Kitaev-k} is highly symmetric. By construction it anticommutes with the particle-hole symmetry $\mathcal{P}=\sigma_x\mathcal{K}$, where $\mathcal{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 $\mathcal{T}=\mathbbm{1}\,\mathcal{K}$. Finally, the product of $\mathcal{TP}=C=\sigma_x$ is the chiral symmetry, whose presence allows us to define the topological invariant in terms of the winding number.~\cite{Review-Chiu-2016} Note that all symmetries square to $+1$, placing the Kitaev chain in the BDI class.~\cite{altland:prb1997}

The winding number
is given by \cite{Wen-1989, Review-Chiu-2016}
\begin{align}\label{equation: winding number invariant}
 	\nu =\frac{1}{2\pi} \int\limits_{-\pi/d}^{\pi/d} dk~\partial_k \,w(k),
\end{align}
where $w(k)=\mathrm{arg}\left[2\Delta\,\sin(kd)\,+\,i\,\left(\mu+2t\,\cos(kd)\right)\right]$ and $\partial_k\, w(k)$ is the winding number density. A trivial phase corresponds to $\nu=0$, a non trivial one to finite integer values of $\nu$. 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 Sec. \ref{section: numerical results and impact of disorder}.

The phase diagram constructed using the winding number invariant is shown in Fig.~\ref{fig:topo-phase-diagram}. 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 \textcolor{blue}{$\gamma_1^A$} connected to the bulk by the $i(\Delta-t)$ hopping and the $B$ Majorana operator \textcolor{orange}{$\gamma_1^B$} connected by the $i(\Delta+t)$ hopping. With $\Delta>0$ and $t<0$ (the $\nu=+1$ phase) the Majorana  state at the left end of the chain will consist mostly of the weakly connected \textcolor{orange}{$\gamma^B_1$}. If $t>0$ (the $\nu=-1$ phase), \textcolor{blue}{$\gamma^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_\pm(k)=0$ (cf. Eq.~\eqref{infinite chain: excitation spectrum}). That is only possible if both terms under the square root vanish. The condition of $\Delta\neq 0$ forces the gap closing to occur at $kd=0$ or $k=\pi d$, and the remaining term vanishes at these momenta if $\mu=\pm2t$. The four insets in Fig. \ref{fig:topo-phase-diagram} show the behavior of $ w (k)$, leading to either a zero (for $w(-\pi)=w(\pi)\,$) or non zero winding number, see Eq. \eqref{equation: winding number invariant}.

Physically speaking, the Kitaev chain is in the topological phase provided that $\Delta\neq 0$ and the chemical potential lies inside the "normal" band ($\vert \mu\vert \le 2 \vert t\vert$).

\section{Spectral analysis of the finite Kitaev chain}
\label{section: spectral analysis}
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 $\Delta$, $t$ and $\mu$ is zero. Already for the simple case $\mu=0$ and $\Delta \neq 0$, $t\neq 0$ the quantization of the momentum turns out to be non trivial. The general case in which all parameters are finite is considered in Secs. \ref{section: influence of non zero chemical potential}, \ref{section: MZM paths in the topological phase diagram for finite N} and \ref{section: numerical results and impact of disorder}.

We start with the BdG Hamiltonian of the open Kitaev chain in real space, expressed in the basis of standard fermionic operators  \makebox{$\hat{\psi}\,=\,\left(d_1,\,\ldots,\,d_N,\,d_1^\dagger,\,\ldots,\,d_N^\dagger\right)^\mathrm{T}$}. Then
\begin{align}\label{equation: Kitaev Hamiltonian and the BdG basis}
	\hat{H}_\mathrm{KC}\,=\,\frac{1}{2}\,\hat{\psi}^\dagger\,\mathcal{H}_\mathrm{KC}\,\psi,
\end{align}
where the BdG Hamiltonian $\mathcal{H}_\mathrm{KC}$ is
\begin{align}\label{equation: BdG Hamiltonian in real space in its full glory}
	\mathcal{H}_\mathrm{KC}\,=\,\left[\begin{matrix}
	C & S\\
	S^\dagger & -C
	\end{matrix}\right].
\end{align}
These matrices have the tridiagonal structure
\begin{align}\label{equation: Matrix C}
	C\,=\,\left[\begin{matrix}
	-\mu & -t  		\\
	-t&-\mu & -t 	\\
	& -t &-\mu & -t &\\
	&  & \ddots & \ddots &\ddots &\\
	&&& -t &-\mu & -t\\
	&&&& -t &-\mu & -t\\
	&&&&& -t &-\mu
	\end{matrix}			\right],
	\\\label{equation: Matrix S}
	S\,=\,\left[\begin{matrix}
	0 & \Delta  		\\
	-\Delta&0 & \Delta 	\\
	& -\Delta &0 & \Delta &\\
	&  & \ddots & \ddots &\ddots &\\
	&&& -\Delta &0 & \Delta\\
	&&&& -\Delta &0 & \Delta\\
	&&&&& -\Delta &0
	\end{matrix}\right].
\end{align}
The spectrum can be obtained by diagonalisation of $\mathcal{H}_\mathrm{KC}$ in real space. We consider different situations.
\subsection{$\Delta=0$}
\label{subsection: spectrum delta  = 0}
The BdG Hamiltonian is block diagonal and its characteristic polynomial $P_\lambda(\mathcal{H}_\mathrm{KC})=\mathrm{det}\left[\lambda\,\mathbbm{1}-\mathcal{H}_\mathrm{KC}\right]$ \makebox{factorises} as
\begin{align}
	P_\lambda(\mathcal{H}_\mathrm{KC})\,=\,P_\lambda(C)\,P_\lambda(-C).
\end{align}
The tridiagonal structure of $C$ straightforwardly yields the spectrum of a normal conducting, linear chain \cite{Kouachi}
\begin{align}\label{spectrum: linear chain}
	E^{\Delta=0}_\pm(k_n)\,=\,\pm\left[\mu\,+\,2t\,\cos\left(k_n d\right)\right],\qquad k_n d\,=\,\frac{n\pi}{N+1},
\end{align}
where $n$ runs from $1$ to $N$. Since $k_n\in\mathbbm{R}$, only bulk states exist for $\Delta = 0$.
\subsection{$t=0$}
\label{subsection: spectrum t  = 0}
In the beginning we consider both $t$ and $\mu$ to be zero and include $\mu\,\neq\,0$ in a second step. The parameter setting leads to a vanishing matrix $C$, see Eq. (\ref{equation: Matrix C}), and the characteristic polynomial of the system reads:
\begin{align}
	P_\lambda(\mathcal{H}_\mathrm{KC})\,=\,\mathrm{det}\left[\begin{matrix}
	\lambda\,\mathbbm{1} & -S\\
	S & \lambda\,\mathbbm{1}
	\end{matrix}\right],
\end{align}
where we used the property $S^\dagger\,=\,-S$. Due to the fact that the commutator $[\mathbbm{1},\,S]=0$ vanishes\footnote{Note that $\lambda$ can be zero and it will for odd $N$. Hence, the standard formula to calculate the determinant of a partitioned $2\times 2$ matrix can not be used here, because it requires the inverse of one diagonal block. We use instead Silvester's formula\cite{Silvester}: $\text{det}\left[\begin{matrix}
A & B\\
C & D\\
\end{matrix} \right] =\text{det} \left[AD\,-\,CB\right]$, where $A,\,B,\,C,\,D$ are square matrices of the same size and the only requirement is $%
\left[C,\,D\right]=0$.}, one finds \cite{Silvester}
\begin{align}
	P_\lambda(\mathcal{H}_\mathrm{KC})\,=\,\mathrm{det}\left(\lambda^2\,\mathbbm{1}\,+\,S^2\right).
\end{align}
The characteristic polynomial can still be simplified to the product
\begin{align}
	P_\lambda(\mathcal{H}_\mathrm{KC})\,=\,P_\lambda(iS)\,P_\lambda(-iS).
\end{align}
The matrix $iS$ is hermitian and describes a {\it{linear}} chain with hopping $i\Delta$. As a consequence, we find the spectrum to be \cite{Kouachi}
\begin{align}\label{equation: spectrum t=mu=0}
	E_\pm(k_n)\,=\,\pm\left[2\Delta\,\cos\left(k_n d\right)\right],\qquad k_n d\,=\,\frac{n\pi}{N+1},
\end{align}
where $n$ runs from $1$ to $N$ and each eigenvalue is twice degenerated. Notice the phase shift by $\pi/2$ compared to the spectrum of an infinite chain Eq. (\ref{infinite chain: excitation spectrum}). We discuss this phase shift in more detail in section \ref{sec:finite-chain_mu0}.

Furthermore if, and only if, $N$ is odd, we find two zero energy modes, namely for \makebox{$n=(N+1)/2$}. Their existence and the degeneracy is due to the chiral symmetry.

The chemical potential $\mu$ can be included easily. Exploiting the properties of $\mathcal{H}_\mathrm{KC}$, we find the characteristic polynomial to be
\begin{align}
	P_\lambda(\mathcal{H}_\mathrm{KC})\,&=\,\mathrm{det}\left[\begin{matrix}
	\left(\lambda +\mu\right)\,\mathbbm{1} & -S\\
	S & \left(\lambda -\mu\right)\,\mathbbm{1}
	\end{matrix}\right]\notag\\
	&=\,\mathrm{det}\left[\left(\lambda^2-\mu^2\right)\,\mathbbm{1}\,+\,S^2\right]\notag\\
	&=\,\mathrm{det}\left[\Lambda^2\,\mathbbm{1}\,+\,S^2\right],
\end{align}
with $\Lambda^2\,\mathrel{\mathop:}=\,\lambda^2-\mu^2$. The same treatment as in the previous $\mu=0$ case yields
$ P_\lambda(\mathcal{H}_\mathrm{KC})\,=\,P_\Lambda(iS)\,P_\Lambda(-iS)$.
Consequently the spectrum is
\begin{align}\label{equation: spectrum and quantization for t  is zero}
	E^{t=0}_\pm(k_n)\,=\,\pm\sqrt{\mu^2\,+\,4\Delta^2\,\cos^2(k_n d)},\qquad k_n d\,=\,\frac{n\pi}{N+1},
\end{align}
where $n$ runs again from $1$ to $N$. Again no boundary modes are found for $t =0$.
\subsection{$\mu=0$}
\label{sec:finite-chain_mu0}
The calculation of the spectrum for $\mu=0$ requires a more technical approach, since the structure of the BdG Hamiltonian Eq. (\ref{equation: BdG Hamiltonian in real space in its full glory}) prohibits standard methods.
\begin{figure}[ht]\centering

	\includegraphics[width = 0.8\columnwidth]{figure2.pdf}
	\caption{Kitaev chain viewed as two coupled SSH-like chains for a) $N=4$ and b) $N=3$ sites. These two chains $\alpha$ and $\beta$ are coupled by $\pm i \mu$. The hoppings $a=i(\Delta-t)$ in red and $b=i(\Delta+t)$ in blue alternate (dashed lines correspond to $-a$ and $-b$) and connect neighbouring Majorana operators $\gamma_j^A$ (blue spheres) and $\gamma_{j\pm 1}^B$ (orange spheres). The unit cell has size $2d$. }
	\label{figure: Kitaev model as coupled ssh chains.}
\end{figure}

One important  feature of the Kitaev chain can be appreciated inspecting Eq. (\ref{equation: Kitaev chain in Majorana operators}). The entire
model is equivalent to two coupled SSH-like chains\cite{wakatsuki:prb2014,li:prb2018} containing both the hopping parameters $a\,\mathrel{\mathop:}=\,i\left(\Delta-t\right)$ and $b\,\mathrel{\mathop:}=\,i\left(\Delta+t\right)$, see Fig. (\ref{figure: Kitaev model as coupled ssh chains.}). Explicitly,
\begin{align}\label{equation: Kitaev Hamilt. in SSH basis, even N}
	\hat{H}_\mathrm{KC}\,&=\,\left(a\sum_{j=1}^{N_1}\,\gamma_{2j-1}^A\gamma_{2j}^B +b \sum_{j=1}^{N_2}\,\gamma_{2j}^B\gamma_{2j+1}^A
	\right)+h.c.\notag\\
	&\,+\left(b\sum_{j=1}^{N_1}\,\gamma_{2j-1}^B\gamma_{2j}^A +a \sum_{j=1}^{N_2}\,\gamma_{2j}^A\gamma_{2j+1}^B \right)+h.c.\notag\\
	&\,-i\mu\sum_{j=1}^N\,\gamma_j^A\gamma_j^B,
\end{align}
where $N_{1,2}$ depend on $N$. If $N$ is even we have $N_1 = N/2$ and $N_2 = N_1-1$, while $N_1=N_2=(N-1)/2$ for odd $N$. Independent of the number of atoms, the first and the second lines in
Eq. (\ref{equation: Kitaev Hamilt. in SSH basis, even N}) describe two SSH-like chains, coupled by the chemical potential $\mu$.
We define here the SSH-like basis of the Kitaev chain as:
\begin{align}\label{equation: ssh basis, alternating pattern}
	\hat{\Psi}^\mathrm{even}_\mathrm{SSH}&=\left(\textcolor{blue}{\gamma_1^\mathrm{A}},\,\textcolor{orange}{\gamma_2^\mathrm{B}},\ldots,					\textcolor{blue}{\gamma_{N-1}^\mathrm{A}},\,\textcolor{orange}{\gamma_N^\mathrm{B}} \vert \textcolor{orange}{\gamma_1^\mathrm{B}},\,				\textcolor{blue}{\gamma_2^\mathrm{A}},\ldots,\textcolor{orange}{\gamma_{N-1}^\mathrm{B}},\,\textcolor{blue}{\gamma_N^\mathrm{A}}	\right)^			\mathrm{T},\notag\\
	\notag\\
	\hat{\Psi}^\mathrm{odd}_\mathrm{SSH}&=\left(\textcolor{blue}{\gamma_1^\mathrm{A}},\,\textcolor{orange}{\gamma_2^\mathrm{B}},\ldots,						\textcolor{orange}{\gamma_{N-1}^\mathrm{B}},\,\textcolor{blue}{\gamma_N^\mathrm{A}} \vert \textcolor{orange}{\gamma_1^\mathrm{B}},\,				\textcolor{blue}{\gamma_2^\mathrm{A}},\ldots,\textcolor{blue}{\gamma_{N-1}^\mathrm{A}},\,\textcolor{orange}{\gamma_N^\mathrm{B}}	\right)^			\mathrm{T},
\end{align}
where "$\vert$" marks the boundary between both chains. We call the first one, starting always with $\gamma_1^\mathrm{A}$, the $\alpha$ chain, and the second the $\beta$ chain, such that \makebox{$\Psi^\mathrm{even,\,odd}_\mathrm{SSH}=\left(\vec{\gamma}_\alpha\,\vert \, \vec{\gamma}_\beta\right)^\mathrm{T}$}. The BdG Hamiltonian in the SSH-like basis reads
\begin{align}\label{equation: Kitaev Hamiltonian/ matrix in SSH basis}
	\mathcal{H}_\mathrm{KC}^\mathrm{SSH}\,=\,\left[\begin{matrix}
	\mathcal{H}_\alpha & \tau \\
	\tau^\dagger &   \mathcal{H}_\beta
	\end{matrix}
	\right],
\end{align}
with $\hat{H}_\mathrm{KC}=\frac{1}{2}\hat{\Psi}_\mathrm{SSH}^\dagger\mathcal{H}_\mathrm{KC}^\mathrm{SSH}\hat{\Psi}_\mathrm{SSH}$. The independent SSH-like chains are represented by the square matrices $\mathcal{H}_\alpha$ and $\mathcal{H}_\beta$ of size $N$. Both chains are coupled by the matrices $\tau$ and $\tau^\dagger$, which contain only the chemical potential $\mu$, 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
\begin{align}
	\mathcal{H}_\alpha^{\mathrm{even}}\,&=\,\left[\begin{matrix}
	0 & a\\
	-a & 0&b \\
	& -b&0& a\\
	&&\ddots & \ddots & \ddots\\
	&&& -a & 0&b\\
	&&&&-b&0& a\\
	&&&&& -a &0
		\end{matrix}		\right],\\
	\mathcal{H}_\beta^{\mathrm{even}}\,&=\,\left[\begin{matrix}
	0 & b\\
	-b & 0&a \\
	& -a&0& b\\
	&&\ddots & \ddots & \ddots\\
	&&& -b & 0&a\\
	&&&&-a&0& b\\
	&&&&& -b &0
		\end{matrix}\,		\right],
\end{align}
and $\tau^{\mathrm{even}} = -i\mu\, \mathbbm{1}_{N/2}\otimes \tau_z$, where $\tau_z$ denotes the Pauli matrix. The odd $N$ expressions are achieved by removing the last line and column in $\mathcal{H}_\alpha^{\mathrm{even}}$, $\mathcal{H}_\beta^{\mathrm{even}}$ and $\tau^\mathrm{even}$.

As shown in more detail in appendix \ref{appendix: spectrum zero mu}, for $\mu=0$ the characteristic polynomial can be expressed as the product of two polynomials of order $N$
\begin{align}\label{equation: characteristic polynomial at zero mu}
	P_\lambda(\mathcal{H}_\mathrm{KC})_{\mu=0}\,=\,\zeta_N(\lambda,\,a,\,b)\,\,\epsilon_N(\lambda,\,a,\,b),
\end{align}
where the product form reflects the fact that the Kitaev chain is given in terms of two uncoupled SSH-like chains, as illustrated in Fig. (\ref{figure: Kitaev model as coupled ssh chains.}). Even though the polynomials $\zeta_N$ and $\epsilon_N$ belong to different SSH-like chains, both obey a common recursion formula typical of Fibonacci polynomials \cite{Webb, Hoggatt, Oezvatan-2017}
\begin{align}\label{equation: recursion formula fibonacci polynomials}
	\zeta_{j+2}\,=\,\left[\lambda^2\,+\,a^2\,+\,b^2\right]\,\zeta_j\,-\,a^2 b^2\,\zeta_{j-2},
\end{align}
and differ only in their initial values
\begin{align}\label{equation: initial values of epsilon and zeta}
	\left(\begin{matrix}
	\zeta_{-1}		\\
	\zeta_{0}		\\
	\zeta_{1}		\\
	\zeta_{2}		\\
\end{matrix}\right)\,=\,\left(\begin{matrix}
	0					\\
	1					\\
	\lambda				\\
	\lambda^2\,+\,b^2
\end{matrix}\right),\,\left(\begin{matrix}
	\epsilon_{-1}		\\
	\epsilon_{0}		\\
	\epsilon_{1}		\\
	\epsilon_{2}		\\
\end{matrix}\right)\,=\,\left(\begin{matrix}
	0					\\
	1					\\
	\lambda				\\
	\lambda^2\,+\,a^2
\end{matrix}\right).
\end{align}
Fundamental properties of Fibonacci polynomials are summarized in appendix \ref{appendix: comments on Fibonacci and Tetrannaci polynomials}. The common sublattice structure of both chains sets the stage for a  relationship between $\zeta_j$ and $\epsilon_j$: The exchange of $a$'s and $b$'s enables us to pass from one to the other
\begin{align}\label{equation: zeta and epsilon, related by abc}
	\zeta_j(\lambda,\,a,\,b)\,=\,\epsilon_j(\lambda,\,b,\,a),\quad \forall j.
\end{align}
Moreover, Eq.\,(\ref{equation: recursion formula fibonacci polynomials}) 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 \cite{Sirker}. The difference emerges since, according to Eqs. (\ref{appendix: closed form zeta_odd}), (\ref{appendix: closed form epsilon_odd}), it holds
\begin{align}\label{equation: odd index polynomials are equal}
	\zeta_\mathrm{odd}(\lambda,\,a,\,b)\,=\,\epsilon_\mathrm{odd}(\lambda,\,a,\,b),
\end{align}
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 $\zeta_j$ and $\epsilon_j$, as well as their factorization, is derived in appendix \ref{appendix: spectrum zero mu}.

The characteristic polynomial can be used to obtain the determinant of the Kitaev chain, here for $\mu=0$, because evaluating it at $\lambda=0$ leads to:
\begin{align*}
	P_{\lambda=0}(\mathcal{H}_\mathrm{KC})_{\mu=0}\,=\,\mathrm{det}\left(\mathcal{H}_\mathrm{KC}\right)_{\mu=0}.
\end{align*}
According to Eq. (\ref{equation: characteristic polynomial at zero mu}) we need only to know $\zeta_N$ and $\epsilon_N$ at $\lambda=0$. The closed form expression for $\zeta_j$ at $\lambda = 0$ reduces to
\begin{align}
	\zeta_j\left\vert_{\lambda=0} \right.\,=\,\left\{\begin{matrix}
	0, &\text{if } j  \text{ is odd}\\
	b^{\,j},&\text{else }
	\end{matrix}\right. ,
\end{align}
while $\epsilon_j \left\vert_{\lambda = 0}\right.$ follows from Eq. (\ref{equation: zeta and epsilon, related by abc}). We find that there are always zero energy eigenvalues for odd $N$, but not in general for even $N$, as it follows from
\begin{align}\label{equation: determinant, pure kitaev chain, N even, mu zero}
	\mathrm{det}\left(\mathcal{H}^{\mu=0}_\mathrm{KC}\right)\,=\,\left\{ \begin{matrix}
	0, & N \text{ odd}\\
	\left[\Delta^2-t^2\right]^N, & N \text{ even}
	\end{matrix}\right. .
\end{align}
Additional features of the spectrum are discussed in the following.
\subsubsection{Odd $N$  }
The spectrum for odd $N$ is given by two contributions
\begin{align}
		\label{equation: spectrum: mu=0, odd N, zero modes}
	E^{\mu=0}_\pm\,&=\,0,\quad \text{(twofold)}\\
		\label{equation: spectrum: mu=0, odd N, k depending part}
	E^{\mu=0}_\pm(k_n )\,&=\,\pm\sqrt{4\Delta^2\,\sin^2(k_n d)\,+\,4t^2\,\cos^2(k_n d)},
\end{align}
where $k_n d=n\pi/(N+1)$ and $n$ runs from $1$ to $N$, except for $n=(N+1)/2$. This constraint on $n$ is a consequence of the Eqs. \eqref{equation: N odd, entries of the eigenstates}, \eqref{equation: function T odd } below which show that the boundary condition Eq. \eqref{equation: boundary condition for odd N ssh chain} cannot be satisfied for $kd=\pi/2$. Hence, no standing wave can be formed. \\
Each zero eigenvalue belongs to one chain. As discussed below, two decaying states are associated to Eq. (\ref{equation: spectrum: mu=0, odd N, zero modes}), whose wave functions are discussed in Sec. \ref{subsection: N odd eigenvectors}. These states are MZM.

\subsubsection{Even $N$}
In the situation of even $N$ we find for the Kitaev's bulk spectrum at zero $\mu$
\begin{align}\label{equation: spectrum: mu=0, even N}
	E^{\mu=0}_\pm(k)\,&=\,\pm\sqrt{4\Delta^2\,\sin^2(k d)\,+\,4t^2\,\cos^2(k d)},
\end{align}
where the momenta $k$ are in general {\it{not equidistant}} in the first Brillouin zone. Rather, the quantization condition follows from the interplay between $\Delta$ and $t$ and is captured in form of the functions $f_{\beta,\alpha}(k)$ (cf. appendix \ref{appendix: factorisation of the fibonnacis}),
\begin{align}\label{Definition f_pm(q)}
	f_{\beta,\alpha}(k)\mathrel{\mathop:}=\tan\left[k d\left(N+1\right)\right]\,\pm\frac{\Delta}{t}\, \tan\left(k d \right) ,
\end{align}
whose zeros
\begin{align}\label{equation: mu=0 quantization condition, N even,}
	f_{\beta,\alpha}(k)\,\stackrel{!}{=}\,0, \quad kd\neq0,\pi/2
\end{align}
define the allowed values of $k$. Note that $kd = 0, \pi/2$ are excluded as solutions, due to their trivial character. The functions $f_{\beta,\alpha}(k)$ follow from the factorisation of the polynomials $\epsilon_N$ and $\zeta_N$.
 The negative sign in Eq. (\ref{Definition f_pm(q)}) belongs to the $\alpha$ subchain, while the positive one to the $\beta$ subchain. The spectrum following from Eq. (\ref{equation: mu=0 quantization condition, N even,}) is illustrated in Fig. \ref{figure: non equidistant bulk quantization for zero mu and N is four}.
\begin{figure}[ht]
	\centering

	\includegraphics[width = 1.0\columnwidth]{figure3.pdf}
	\caption{Eigenvalues and the non equidistant quantization of the bulk momentum $k$ for a Kitaev molecule with four sites. a) The horizontal lines mark the numerical eigenvalues $\pm E_{j}$ ($j=0,\,1,\,2,\,3$) and the bulk spectrum of the infinite chain is the green solid
 curve. The tangent-like functions follow $f_\beta(k)$ and $f_\alpha(k)$ from Eq. (\ref{Definition f_pm(q)}). The situation shown in a) is for $k\in\left[-\pi/d, \,\pi/d\right]$. The zeros of $f_{\beta,\alpha}(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_\pm(k_j)=\pm E_{j}$. b) Zoom of a) for $k\in\left[0.5/d,\,1.2/d\right]$. The chosen parameters $N=4$, $t=4\,$meV, $\Delta = 1.5\,$ meV and $\mu = 0$ meV lead to the bulk eigenvalues $\pm E_{j}\in\left[\, \pm 4.39,\,\pm 6.47,\,\pm 6.89\right]$ (in meV) and to the momenta $k_{3,\,2,\,1}$ approximately $[0.58012/d, 0.68813/d, 1.12386/d]$.}
	 \label{figure: non equidistant bulk quantization for zero mu and N is four}
\end{figure}
We observe that Eqs. (\ref{equation: spectrum: mu=0, even N}) and (\ref{equation: mu=0 quantization condition, N even,}) hold for all values of $t$ and $\Delta$, independent of whether $\vert \Delta\vert$ is larger or smaller than $\vert t \vert$. The two situations are connected by a phase shift of the momentum $kd\rightarrow kd+\pi/2$, which influences both the spectrum and the quantization condition. In the end all different ratios of $\Delta$ and $t$ are captured by Eqs. (\ref{equation: spectrum: mu=0, even N}) and (\ref{equation: mu=0 quantization condition, N even,}), due to the periodicity of the spectrum.

However, when we consider decaying or edge states this periodicity is lost (see Eqs. (\ref{equation: criteria/ number of solutions for edge states even N}) - (\ref{equation: spectrum: mu=0, even N, edges}) below) and $\vert t\vert\lessgtr \vert \Delta\vert$ lead to different quantization rules.
The hermiticity of the Hamiltonian allows a pure imaginary momentum for $\mu=0$, but a simple exchange of $k$ to $i q$ in Eq. \eqref{Definition f_pm(q)} does not lead to the correct results. We introduce here the functions
\begin{align}\label{definiton: h_pm(q)}
	h_{\beta,\alpha}(q)\mathrel{\mathop:}=\tanh\left[q d\left(N+1\right)\right]\,\pm\,m\,\tanh\left(q d\right),
\end{align}
similar to $f_{\beta,\alpha}(k)$ in Eq. (\ref{Definition f_pm(q)}), where $m$ contains both ratios of $\Delta$ and $t$:
\begin{align}\label{equation: criteria of the prop. factor m}
	m\,\mathrel{\mathop:}= \,\left\{	\begin{matrix}
	\frac{\Delta}{t},\quad\,\mathrm{if}\, \vert \Delta\vert \ge \vert t\vert	\\
	\\
	\frac{t}{\Delta},\quad\,\mathrm{if}\, \vert t\vert \ge \vert \Delta\vert
	\end{matrix}		\right. .
\end{align}
Again, the positive sign in Eq. (\ref{definiton: h_pm(q)}) belongs to the $\beta$ chain and the negative one to the $\alpha$ chain. The exact quantization criterion is provided by the zeros of $h_{\beta,\alpha}(q)$,
\begin{align}\label{equation: criteria/ number of solutions for edge states even N}
	h_{\beta,\alpha}(q)\,\stackrel{!}{=}0, \quad q\neq 0,
\end{align}
as illustrated in Fig. \ref{figure: non equidistant decaying states quantization for zero mu and N is four}. The associated energies follow from the dispersion relation
\begin{align}\label{equation: spectrum: mu=0, even N, edges}
	E(q)=\pm\,\left\{\begin{matrix}
	\sqrt{4t^2\cosh^2(qd)-4\Delta^2\,\sinh^2(qd)}, \,\mathrm{if}\, \vert \Delta\vert \ge \vert t\vert	\\
	\\
	\sqrt{4\Delta^2\cosh^2(qd)-4t^2\,\sinh^2(qd)},\,\mathrm{if}\, \vert t	\vert \ge \vert \Delta \vert
	\end{matrix}			\right..
\end{align}
We notice that Eq. (\ref{equation: spectrum: mu=0, even N, edges}) is only well defined for zero or positive arguments  of the square root. Indeed, all solutions of Eq. (\ref{equation: criteria/ number of solutions for edge states even N}), if existent, lie always inside this range, because using $h_{\beta,\alpha}(q) = 0$ in Eq. (\ref{equation: spectrum: mu=0, even N, edges}) yields
\begin{align}\label{equation: spectrum: mu=0, even N, edges, cosh/cosh version}
	E(q)\,=\,\pm\,2\,\frac{\cosh(qd)}{\cosh\left[qd\left(N+1\right)\right]}\,\cdot\,\mathrm{min}\left\{\vert \Delta\vert,~ \vert t \vert\right\}.
\end{align}
Hence, each wavevector from Eq. (\ref{equation: criteria/ number of solutions for edge states even N}) \makebox{corresponds} to two gap modes, since the gap width is $4\,\mathrm{min}\left\{\vert \Delta\vert,~ \vert t \vert\right\}$ and the fraction inside Eq. (\ref{equation: spectrum: mu=0, even N, edges, cosh/cosh version}) is always smaller than one.
\begin{figure}[ht]
	\centering

	\includegraphics[width = 1.0\columnwidth]{figure4.pdf}
	\caption{Eigenvalues and the quantised momentum $q_0$ of the gap modes for a Kitaev molecule with four sites. a) The horizontal lines characterise the numerical eigenvalues $\pm E_{j}$ ($j=0,\,1,\,2,\,3$) and the dispersion relation inside the gap is shown as function of continuous $q$ on a finite range. Only one of both hyperbolic tangent-like functions $h_\beta(q)$ or $h_\alpha(q)$ defines a proper $qd\neq 0$. The situation shown in a) is for $q\in\left[-\pi/3d, \,\pi/3d\right]$. b) Zoom of a) for $q\in\left[0.3/d,\,0.42/d\right]$. The momenta $\pm q_0$, the zero of $h_\alpha$, leads to the correct associated energies, such that $E_{\pm}(q_0)=\pm E_{0}$. The chosen parameters $N =4$, $t=4\,$meV, $\Delta = 1.5\,$ meV and $\mu = 0$ meV lead to the eigenvalues $E_{j,\pm}\in\left[\pm 0.97,\, \pm 4.39,\,\pm 6.47,\,\pm 6.89\right]$ (in meV) and to the momentum $q_0\approx 0.37416/d $.}
	 \label{figure: non equidistant decaying states quantization for zero mu and N is four}
\end{figure}
We can restrict ourselves to find only positive solutions $qd$, due to the time reversal symmetry. The number of physically different solutions
of Eq. (\ref{equation: criteria/ number of solutions for edge states even N}) is zero or two and it follows always from the equation containing the positive factor $m$ or $-m$. Consequently, according to Eq. (\ref{definiton: h_pm(q)}), only zero or two gap modes can form and both belong to the same subchain, $\alpha$ or $\beta$. Moreover a solution exists if, and only if, $\vert m\vert \in\left[1,\,N+1\right]$.

In the limiting case when \makebox{$\vert m \vert\rightarrow 1$}, i.e. at the Kitaev points, the solution \makebox{$qd \rightarrow \infty$} and the associated energies $E_\pm$ from Eq. (\ref{equation: spectrum: mu=0, even N, edges, cosh/cosh version}) go to zero. The eigenstate will be a Majorana zero energy mode, see Sec. \ref{subsection: N even eigenvectors}.

In the second special case of $\vert m \vert\rightarrow N+1$ the solution approaches zero. The value $q=0$ is only in this particular scenario a proper momentum, see appendix \ref{appendix: factorisation of the fibonnacis}.
The momentum $q = 0$ yields the energies $E_\pm(0) = \pm 2\,\mathrm{min}\left\{\vert \Delta\vert,~ \vert t \vert\right\}$, which mark exactly the gap boundaries.

Increasing the value of $\vert m \vert$ 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 Eq. (\ref{equation: mu=0 quantization condition, N even,}) 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 $\xi\propto 1/q$ reached the highest possible delocalisation in the chain.

The limit of $N\rightarrow\infty$ yields always two zero energy boundary modes; since the momentum is $qd=\text{arctanh}(\vert 1/ m \vert)$, due to Eqs. (\ref{definiton: h_pm(q)}) (\ref{equation: criteria/ number of solutions for edge states even N}) and according to Eq. (\ref{equation: spectrum: mu=0, even N, edges, cosh/cosh version}) 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 \ref{table: all quantisiation (sub) rules}.
\section{Eigenvectors ($\mu = 0$)}
\label{section: Eigenvectors and symmetries}
We use the SSH-like basis to calculate the eigenvectors of the Hamiltonian Eq. (\ref{equation: Kitaev Hamiltonian/ matrix in SSH basis}) at $\mu=0$. The eigenvectors $\vec{\psi}$ are defined with respect to the SSH-like chains $\alpha$ and $\beta$, see Eq. (\ref{equation: Kitaev Hamiltonian/ matrix in SSH basis}),
\begin{align}
	\vec{\psi}\,=\,\left(\begin{matrix}
	\vec{v}_\alpha\\
	\vec{v}_\beta
\end{matrix}\right),
\end{align}
with the feature that always either $\vec{v}_\beta$ or $\vec{v}_\alpha$ can be chosen to be zero, yielding the solutions $\vec{\psi}^\alpha$ and $\vec{\psi}^\beta$, respectively
\begin{align}\label{equation: eigenvector structure psi_alpha, psi_beta}
\vec{\psi}^\alpha \,=\,\left(\begin{matrix}
	\vec{v}_\alpha\\
	\vec{0}
\end{matrix}\right),\quad\vec{\psi}^{\beta} \,=\,\left(\begin{matrix}
	\vec{0}\\
	\vec{v}_\beta
\end{matrix}\right).
\end{align}
We are left to find the eigenvectors of a single tridiagonal matrix which we did basing on, and extending the results of Ref.~[\onlinecite{Shin-97}]. We focus here on the edge and decaying states, while the rest of our results are in appendix \ref{appendix: eigenvectors}. Remember that in the SSH-like basis Eq. (\ref{equation: ssh basis, alternating pattern}) the Majorana operators $\gamma_j^A$ and $\gamma_j^B$, alternate at each site, thus defining two interpenetrating "A" and "B" type sublattices.

\subsection{Even $N$}
\label{subsection: N even eigenvectors}
We define the vectors $\vec{v}_\alpha$ and $\vec{v}_\beta$ via the entries
\begin{align}
	\vec{v}_\alpha\,&=\,\left(x_1,\,y_1,\,x_2,\,y_2,\,\ldots ,\,x_{N/2},\,y_{N/2}\right)^\mathrm{T},\\
	\vec{v}_\beta\,&=\,\left({\scriptstyle{\mathcal{X}}}_1,\,{\scriptstyle{\mathcal{Y}}}_1,\,{\scriptstyle{\mathcal{X}}}_2,\,{\scriptstyle{\mathcal{Y}}}_2,\,\ldots ,\,{\scriptstyle{\mathcal{X}}}_{N/2},\,{\scriptstyle{\mathcal{Y}}}_{N/2}\right)^\mathrm{T},
\end{align}
where $x,\,y$ and ${\scriptstyle{\mathcal{X}}},\,{\scriptstyle{\mathcal{Y}}}$ are associated to the A and B sublattices, respectively. The internal structure of $\vec{v}_\alpha$ ($\vec{v}_\beta$) reflects the unit cell of an SSH-like chain and thus simplifies the calculation.
\begin{figure}[ht]

	\includegraphics[width = 1\columnwidth]{figure5.pdf}
	\caption{Visualisation of the entries of the eigenstates $\vec{\psi}^\alpha_{q,k} (j)$ of the Kitaev chain with $N =42$ sites and $\mu =0$. Panel a) depicts the gap state $\vec{\psi}^\alpha_q$ and b) the lowest energy bulk state $\vec{\psi}^\alpha_k$. The blue (orange) dots follow $x_l/x_1$ ($i\,y_l/x_1$) at position $j = 2l-1$ ($j= 2l$) for $l=1,\ldots, N/2$, while the black line is only a guide to the eye. The gap state is more localised at the edges. The extended state is largest inside the chain. The chosen parameters are $t=10$ meV and $\Delta = 1$ meV leading to $q = 0.10029/d$ and $E = 0.0539$ meV for the gap state. The shown extended state is associated with $k =1.4806/d$ and $E = 2.6851$ meV. Notice the decaying state in a) as well as the ones depicted in Figs. \ref{figure: wavefunctions N even, t=5} a) and \ref{figure: wavefunctions N even, t=-5, beta chain} are not Majorana states.
	}
	\label{figure: wavefunctions N even, t=10}
\end{figure}
In the real space  $x_l$ (${\scriptstyle{\mathcal{X}}}_l$) belongs to site $j = 2l-1$ and $y_l$ (${\scriptstyle{\mathcal{Y}}}_l$) to $j=2l$, where $j=1,\,\ldots,\,N$.
\begin{figure}[h]

	\includegraphics[width = 1\columnwidth]{figure6.pdf}
	\caption{Illustration of the entries of the eigenstates $\vec{\psi}^\alpha_{q,k}(j) $ of a Kitaev chain for $N=42$ sites and $\mu=0$ for $t =5$ meV and $\Delta = 1$ meV. Similar to Fig. \ref{figure: wavefunctions N even, t=10}, but with a modified value for $t/\Delta$. a) shows the gap mode and b) the lowest in energy bulk mode. Notice that for the chosen parameter set the gap state is more localized than the one in Fig. \ref{figure: wavefunctions N even, t=10}. 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 $E =0.6682\cdot 10^{-3}$ meV ($E= 2.1555$ meV). }
	\label{figure: wavefunctions N even, t=5}
\end{figure}

Searching for solutions on the subchain $\alpha$ implies setting \makebox{$\vec{v}_\beta = \vec{0}$} and solving $\left(\mathcal{H}_\alpha^{\mathrm{even}  } -E_\pm \mathbbm{1}_N\right)\vec{v}_\alpha = \vec{0}$. The elements of $\vec{v}_\alpha$ obey
\begin{align}
	 \label{equation: open boundary left side}
	 \,a\,y_1\,&=\,E_\pm\,x_1,\\
	 \label{equation: open boundary right side}
	-a\,x_{N/2}\,&=\,E_\pm\,y_{N/2},
\end{align}
and
\begin{align}
	 \label{equation: recursive formula, entries v_alpha x}
	 b\,x_{l+1}\,-\,a\,x_l\,&=\,E_\pm\,y_l,\\
	 \label{equation: recursive formula, entries v_alpha y}
	 a\,y_{l+1}\,-\,b\,y_l\,&=\,E_\pm\,x_{l+1},
\end{align}
where $l$ runs from $1$ to $N/2 \,-\, 1$. The solution for $\Delta \neq \pm t$ is (in agreement with Ref.~[\onlinecite{Shin-97})]
\begin{align}
	\label{equation: decaying state, even N, y_j}
	\frac{y_l}{x_1}\,&=\,\frac{E_\pm (\theta)}{a}\,T^\mathrm{e}_l(\theta),\\
	\label{equation: decaying state, even N, x_j}
	\frac{x_l}{x_1}\,&=\,T^\mathrm{e}_l(\theta)-\frac{b}{a} T^\mathrm{e}_{l-1}(\theta),
\end{align}
where $l = 1,\,\ldots,\,N/2$, $\theta$ denotes the momentum $k$ ($q$) for extended (gap) states and $E_\pm$ is the dispersion relation associated to $k$ (Eq. (\ref{equation: spectrum: mu=0, even N})), or $q$ (Eq. (\ref{equation: spectrum: mu=0, even N, edges})).
\begin{figure}

		\includegraphics[width = 1\columnwidth]{figure7.pdf}
	\caption{The decaying state $\vec{\psi}^\beta_q$ for $N=42$ sites and $\mu=0$. The black guiding line follows the orange (blue) dots, which correspond to ${\scriptstyle{\mathcal{X}}}_l/{\scriptstyle{\mathcal{X}}}_1$  ($i\,{\scriptstyle{\mathcal{Y}}}_l/{\scriptstyle{\mathcal{X}}}_1$) at position $j=2l-1$ ($j= 2l$), $l = 1,\ldots, N/2$. The difference to the edge states on subchain $\alpha$ is the exchanged role of Majorana operators $\gamma_j^A$, $\gamma_j^B$. The chosen parameters are $t=-5$ meV and $\Delta = 1$ meV. The gap mode is associated with $q = 0.2027/d$ and $E =0.6682\cdot 10^{-3}$ meV .  }
	\label{figure: wavefunctions N even, t=-5, beta chain}
\end{figure}
The entries of the eigenvectors are essentially sine functions for the extended states
\begin{align}\label{equation: T even, bulk states}
	T^\mathrm{e}_l (k)\mathrel{\mathop:}=\,\frac{\sin(2\,kd\,l)}{\sin(2\,kd)},
\end{align}
and hyperbolic sine functions for the decaying states
\begin{align}\label{equation: T even, gap modes}
	T^\mathrm{e}_l (q)\mathrel{\mathop:}= s^{l-1}\,\,\,\frac{\sinh(2\,qd\,l)}{\sinh(2\,qd)},
\end{align}
where the prefactor $s$ depends on the ratio of $\Delta$ and $t$:
\begin{align*}
	s = \left\{\begin{matrix}
	+1, & \vert \Delta\vert > \vert t\vert\\
	\\
	-1, & \vert t\vert > \vert \Delta\vert
	\end{matrix}\right. .
\end{align*}
An illustration of $\vec{\psi}^\alpha$ is given in Fig. \ref{figure: wavefunctions N even, t=10}. The allowed momenta $k$ or $q$ follow from the open boundary conditions
\begin{align}
	\label{equation: boundary condition y_0, x_N/2+1 has to vanish}
	y_0\,=\,x_{\frac{N}{2}+1}\,=\,0.
\end{align}
The first condition is satisfied due to $T_0(\theta)=0$ for any momentum. The second condition yields the quantization rules $f_\alpha(k)=0$ and $h_\alpha(q)=0$ for the $\alpha$ chain, see Eqs. (\ref{equation: mu=0 quantization condition, N even,}), (\ref{equation: criteria/ number of solutions for edge states even N}).

The eigenvector $\vec{\psi}^\beta$ entails $\vec{v}_\alpha=0$ and the entries of $\vec{v}_\beta$ follow essentially by replacing $a$'s
and $b$'s in the Eqs. (\ref{equation: decaying state, even N, y_j}), (\ref{equation: decaying state, even N, x_j}). We find
\begin{align}\label{equation: decaying state on chain beta, even N, ycal_j}
	\frac{{\scriptstyle{\mathcal{Y}}}_l}{{\scriptstyle{\mathcal{X}}}_1}\,&=\,\frac{E_\pm }{b}\,T^\mathrm{e}_l(\theta),\\
	\label{equation: decaying state on chain beta, even N, xcal_j}
	\frac{{\scriptstyle{\mathcal{X}}}_l}{{\scriptstyle{\mathcal{X}}}_1}\,&=\,T^\mathrm{e}_l(\theta)-\frac{a}{b} T^\mathrm{e}_{l-1}(\theta),
\end{align}
where $j = 1,\,\ldots,\,N/2$ and $\Delta \neq \pm t$. The quantisation condition follows from the open boundary condition:
\begin{align*}
	{\scriptstyle{\mathcal{Y}}}_0\,&=\,0,\quad {\scriptstyle{\mathcal{X}}}_{\frac{N}{2}+1}\,=\,0,
\end{align*}
and $k$ ($q$) obey $f_\beta(k)=0$ ($h_\beta(q)=0$). Further, from the quantization rules it follows that gap modes belong always to the same subchain $\alpha$ or $\beta$ for even $N$.

As illustrated in Figs. \ref{figure: wavefunctions N even, t=10}, \ref{figure: wavefunctions N even, t=5}, \ref{figure: wavefunctions N even, t=-5, beta chain} our states are symmetric w.r.t. the center of the SSH-like chains. This symmetry is visible in alternative versions of the Eqs. \eqref{equation: decaying state, even N, x_j}, \eqref{equation: decaying state on chain beta, even N, xcal_j} ($l=1,\ldots,N/2$), whereby
$	{\scriptstyle{\mathcal{X}}}_{\frac{N}{2}+1-l}\,=\,T_l^\mathrm{e}(\theta)\,{\scriptstyle{\mathcal{X}}}_{\frac{N}{2}}$,
 which holds for all eigenstates. Together with Eqs. \eqref{equation: decaying state, even N, y_j},  we find in general
\begin{align}\label{equation: inversion symmetry of the even N ssh like chain for mu=0}
	x_{\frac{N}{2}+1-l}=y_{l}\,\frac{x_{N/2}}{y_1}
\end{align}
and similarly ${\scriptstyle{\mathcal{X}}}_{\frac{N}{2}+1-l}={\scriptstyle{\mathcal{Y}}}_{l}\,{\scriptstyle{\mathcal{X}}}_{\frac{N}{2}}/{\scriptstyle{\mathcal{Y}}}_1$. Recalling the definition of the SSH-like basis, Eq. \eqref{equation: ssh basis, alternating pattern}, and introducing the operators $\psi_{\alpha,\beta}^\dagger$ associated to the states $\vec{\psi}^{\alpha,\beta}$ in Eq. \eqref{equation: eigenvector structure psi_alpha, psi_beta}, we find the expression
\begin{align}\label{equation: operator alpha N even}
	\psi_{\alpha}^\dagger\,=\,\frac{1}{v_\alpha}\left[\sum\limits_{j=1}^{N/2}\,x_j\,\left(\gamma_{2j-1}^A\right)^\dagger\,+\,\sum\limits_{j=1}^{N/2}\,y_j\,\left(\gamma_{2j}^B\right)^\dagger\right],
\end{align}
where $v_\alpha$ is the norm of the vector $\vec{v}_\alpha$. A similar term is found for $\psi_{\beta}^\dagger$.
We notice that Eq. \eqref{equation: operator alpha N even} and Eq. \eqref{equation: statistics of the even N eigenstates} below are true for all kinds of eigenstates, i.e. extended, decaying states and MZM, of the BdG Hamiltonian in Eq. \eqref{equation: Kitaev Hamiltonian/ matrix in SSH basis} at $\mu=0$. The character (statistics) of the operators depends on whether $(\psi_{\alpha,\beta}^\dagger)^2$ is $0$ or $1/2$. The property $\{\gamma^r_j,\gamma^s_k\}=\delta_{j,k}\delta_{r,s}$ with ($r,s \in \{A,B\}$) yields
\begin{align}\label{equation: statistics of the even N eigenstates}
	\left(\psi_{\alpha}^\dagger\right)^2\,=\,\frac{1}{2\,v_\alpha^2}\sum\limits_{j=1}^{N/2}\,\left(x_j^2\,+\,y_j^2\right).
\end{align}
The symmetry in Eqs. \eqref{equation: inversion symmetry of the even N ssh like chain for mu=0} states that $(\psi_{\alpha}^\dagger)^2$ is essentially determined by $(x_{N/2}/y_1)^2$. For $a\neq0$ and consequently $E\neq 0$, we find from Eqs. \eqref{equation: open boundary left side}, \eqref{equation: open boundary right side} and \eqref{equation: inversion symmetry of the even N ssh like chain for mu=0} that $(x_{N/2}/y_1)^2 = -1$, which yields $(\psi_{\alpha}^\dagger)^2 =0$. Thus the operators associated to the finite energy states $\vec{\psi}^\alpha$, including the ones depicted in Figs. \ref{figure: wavefunctions N even, t=10} a) and \ref{figure: wavefunctions N even, t=5} a), obey fermionic statistics. This result holds also true in the case $a=0$ and $E\neq 0$ as can be seen by using the corresponding eigenstates (appendix \ref{appendix: eigenvectors}). Similar results hold for $(\psi_{\beta}^\dagger)^2$.

We turn now to Majorana zero modes, which at $\mu=0$ only exist at the Kitaev points $\Delta = \pm t$.

When $\Delta = t$ we find two zero energy modes $\vec{\psi}_A^\alpha = \left(\begin{matrix}
\vec{v}_{\alpha, A}\\ \vec{0}\end{matrix}\right)$, $\vec{\psi}_B^\alpha = \left(\begin{matrix}
\vec{v}_{\alpha, B}\\ \vec{0}\end{matrix}\right)$ each localised at one end of the $\alpha$ chain:
\begin{align}\label{equation: Delta = t, N even}
	\vec{v}_{\alpha,\,\textcolor{blue}{A}}\,&=\,\left(1,\,0,\,0,\,\ldots,\,0\right)^\mathrm{T},\\
	\vec{v}_{\alpha,\,\textcolor{orange}{B}}\,&=\,\left(0,\,0,\,\ldots,\,0,\,1\right)^\mathrm{T},
\end{align}
and $(\psi_\alpha^\dagger)^2 = 1/2$ in Eq. \eqref{equation: statistics of the even N eigenstates}. In contrast, both zero energy modes are on the $\beta$ chain for \makebox{$\Delta = -t$}. We find $\vec{\psi}_A^\beta = \left(\begin{matrix} \vec{0}\\
\vec{v}_{\beta, A} \end{matrix}\right)$, $\vec{\psi}_B^\beta = \left(\begin{matrix}\vec{0}\\
\vec{v}_{\beta, B} \end{matrix}\right)$ with
\begin{align}
	\vec{v}_{\beta,\,\textcolor{orange}{B}}\,&=\,\left(1,\,0,\,0,\,\ldots,\,0\right)^\mathrm{T},\\
	\label{equation: Delta = -t, N even}
	\vec{v}_{\beta,\,\textcolor{blue}{A}}\,&=\,\left(0,\,0,\,\ldots,\,0,\,1\right)^\mathrm{T}.
\end{align}
These states are the archetypal Majorana zero modes\cite{kitaev:physusp2001, Aguado}. Due to their degeneracy, these modes can be recombined into fermionic quasiparticles by appropriate linear combination, see in Eq. \eqref{equation: Missing electron a}, \eqref{equation: Missing electron b} from Sec. \ref{section 2 model}.
\subsection{Odd $N$ }
\label{subsection: N odd eigenvectors}
The composition of the eigenvectors slightly changes for the odd case compared to the even $N$ case
\begin{align}
	\vec{v}_\alpha\,&=\,\left(x_1,\,y_1,\,x_2,\,y_2,\,\ldots ,\,x_{\frac{N-1}{2}},\,y_{\frac{N-1}{2}},\,x_{\frac{N+1}{2}}\right)^\mathrm{T},\\
	\vec{v}_\beta\,&=\,\left({\scriptstyle{\mathcal{X}}}_1,\,{\scriptstyle{\mathcal{Y}}}_1,\,{\scriptstyle{\mathcal{X}}}_2,\,{\scriptstyle{\mathcal{Y}}}_2,\,\ldots ,\,{\scriptstyle{\mathcal{X}}}_{\frac{N-1}{2}},\,{\scriptstyle{\mathcal{Y}}}_{\frac{N-1}{2}},\,{\scriptstyle{\mathcal{X}}}_{\frac{N+1}{2}}\right)^\mathrm{T}.
\end{align}
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 $(\Delta\,\neq\,\pm t$ \emph{and} $E_\pm \neq 0$) does not differ much from the one of the even chain and the entries of $\vec{v}_\alpha$ are
\begin{align}\label{equation: N odd, entries of the eigenstates}
	\frac{y_l}{x_1}\,&=\,\frac{E_\pm(k_n)}{a}\,T^\mathrm{o}_l(k_n ),\\
	\label{equation: N odd, entries of the eigenstates, x_l}
	\frac{x_l}{x_1}\,&=\,T^\mathrm{o}_l(k_n )\,-\,\frac{b}{a}\,T^\mathrm{o}_{l-1}(k_n ),
\end{align}
where $T^\mathrm{o}_l$ is
\begin{align}\label{equation: function T odd }
	T^\mathrm{o}_l (k_n )\mathrel{\mathop:}= \frac{\sin(2\,k_n d\,l)}{\sin(2\,k_n d)},
\end{align}
with $k_n d\,=\,n\pi/(N+1)$ ($n=1,\,\ldots,\,N$, $n\neq (N+1)/2$). The exchange of $a$'s and $b$'s leads again to the coefficients for the chain $\beta$ (see appendix \ref{appendix: eigenvectors}).

The significant difference between even and odd $N$ lies in the realization of the open boundary condition. Solving $\left(\mathcal{H}_\alpha^{\mathrm{odd}  } -E_\pm \mathbbm{1}_N\right)\vec{v}_\alpha = \vec{0}$ yields now
\begin{align}\label{equation: boundary condition for odd N ssh chain}
	y_0\,=\,0,\quad 	y_{\frac{N+1}{2}}\,=\,0,
\end{align}
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 $\alpha$ and $\beta$ contribute each with one.
We find on subchain $\alpha$ for $\Delta\neq\pm t$
\begin{align}\label{equation: zero mode on alpha N odd}
	y_l\,=\,0,\quad 	x_l\,=\,\left(\frac{\Delta-t}{\Delta+t}\right)^{l-1}\,x_1,
\end{align}
and on subchain $\beta$
\begin{align}\label{equation: zero mode on beta N odd}
	{\scriptstyle{\mathcal{Y}}}_l\,=\,0,\quad
	{\scriptstyle{\mathcal{X}}}_l\,=\,\left(\frac{\Delta+t}{\Delta-t}\right)^{l-1}\,{\scriptstyle{\mathcal{X}}}_1,
\end{align}
where $l$ runs from $1$ to $(N+1)/2$.

Regarding the statistics of the operators $\psi_{\alpha}^\dagger,\,\psi_{ \beta}^\dagger$ associated to the states $\vec{\psi}^{\alpha},\, \vec{\psi}^{\beta}$, we proceed like for the even $N$ case. The use of the SSH-like basis from Eq. \eqref{equation: ssh basis, alternating pattern} and the entries of the state $\vec{\psi}^\alpha$ yield now
\begin{align*}
	\psi_{\alpha}^\dagger\,=\,\frac{1}{v_\alpha}\left[\sum\limits_{j=1}^{\frac{N+1}{2}}\,x_j\,\left(\gamma_{2j-1}^A\right)^\dagger\,+\,\sum\limits_{j=1}^{\frac{N-1}{2}}\,y_j\,\left(\gamma_{2j}^B\right)^\dagger\right].
\end{align*}
Again, the Eqs. \eqref{equation: N odd, entries of the eigenstates}, \eqref{equation: N odd, entries of the eigenstates, x_l} and  \eqref{equation: boundary condition for odd N ssh chain} imply a perfect compensation of the $A$ and $B$ sublattice contributions, yielding $(\psi_{\alpha}^\dagger)^2=0$ for $E\ne0$. The zero energy mode, given by its entries in Eq. \eqref{equation: zero mode on alpha N odd}, leads to $(\psi_{\alpha}^\dagger)^2=1/2$.

Further, we find that both zero energy modes $\vec{\psi}^{\alpha, \beta}$ 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
\begin{align}\label{equation: decay length}
	\xi= 2d \left\{ \begin{array}{cc}
	                \left\vert \ln\left(\frac{t-\Delta}{t+\Delta}\right)\right\vert^{-1} & \vert t\vert \geq \vert \Delta \vert, \\[2mm]
	                \left\vert \ln\left(\frac{\Delta-t}{\Delta+t}\right)\right\vert^{-1} & \vert t\vert \leq \vert \Delta \vert,
	                \end{array}\right.
	                \end{align}
and remembering that the atomic site index of $x_l$ is $j=2l-1$ Eq. (\ref{equation: zero mode on alpha N odd}) yields for $\vert t\vert \geq \vert \Delta \vert$
\begin{align}
	x_l  \,&=\, x_1\, (-1)^{l-1}\, e^{-2(l-1)d/\xi},\notag\\
	&=\,x_1\, (-1)^{l-1}\, e^{-(j-1)d/\xi}.
\end{align}
For $\vert t\vert \leq \vert \Delta \vert$ the $x_l$ coefficients are given by the same equation without the $(-1)^{l-1}$ factor. We have moreover $q=\pm1/\xi$, where $q$ is the imaginary momentum yielding $E=0$ in Eq.~\eqref{equation: spectrum: mu=0, even N, edges}.
Thus the localisation of these states is determined only by $t$ and $\Delta$. In the parameter setting of $\Delta =t$ we find:
\begin{align}
	\label{equation: zero mode on alpha, odd N, Kitaev point plus}
	\vec{v}_{\alpha,\,\textcolor{blue}{A}}\,&=\,\left(1,\,0,\,0,\,\ldots,\,0\right)^\mathrm{T},\\
	\vec{v}_{\beta,\,\textcolor{orange}{B}}\,&=\,\left(0,\,0,\,\ldots,\,0,\,1\right)^\mathrm{T},
\end{align}
while both states exchange their position for \makebox{$\Delta = -t$}
\begin{align}
	\vec{v}_{\alpha,\,\textcolor{orange}{B}}\,&=\,\left(1,\,0,\,0,\,\ldots,\,0\right)^\mathrm{T},\\
	\label{equation: zero mode on beta, odd N, Kitaev point minus}
	\vec{v}_{\beta,\,\textcolor{blue}{A}}\,&=\,\left(0,\,0,\,\ldots,\,0,\,1\right)^\mathrm{T}.
\end{align}
\subsection{The particle-hole-operator}
\label{subsection: three symmetries}
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 $\psi_{\alpha,\beta}$. We further corroborate this statement now by recalling that a MZM is defined as an eigenstate of the Hamiltonian $\mathcal{H}$ and of the particle hole symmetry $\mathcal{P}$. The latter acts on an eigenstate $\vec{\psi}^{\alpha, \beta}$ of energy $E$ by turning it into an eigenstate of $\mathcal{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 particle-hole symmetry, discussed in Sec. \ref{subsection: Topological phase diagram}, can be constructed in real space too. Of particular interest is their representation in the SSH-like basis. The antiunitary particle-hole symmetry is
\begin{align}\label{equation: particle hole symmetry in real space in ssh basis}
	\mathcal{P}\,=\,\mathcal{K}\,\mathbbm{1}.
\end{align}
The time reversal and the chiral symmetry depend on $N$. If $N$ is even we find
\begin{align}
	\label{equation: chiral symmetry even N}
	\mathcal{C}^{\,\mathrm{even}}\,&=\,\left[\begin{matrix}
	\mathbbm{1}_{N/2}\otimes \tau_z & \\
	& -\mathbbm{1}_{N/2}\otimes \tau_z
	\end{matrix}
	\right],\\
	\mathcal{T}^{\,\mathrm{even}}\,&=\,\left[\begin{matrix}
	\mathbbm{1}_{N/2}\otimes \tau_z & \\
	& -\mathbbm{1}_{N/2}\otimes \tau_z
	\end{matrix}
	\right]\,\mathcal{K}.
\end{align}
The expressions for odd $N$ follow by removing the last line and last column in each diagonal block.

The effect of $\mathcal{P}$ from Eq. (\ref{equation: particle hole symmetry in real space in ssh basis}) can be seen explicitly if one considers $\mathcal{P}\,\vec{\psi}^{\alpha}$. For $N$ even and $\Delta \neq t$
 the elements $x_l$, $y_l$ of $\vec{\psi}^{\alpha}$ are given in Eqs. (\ref{equation: decaying state, even N, y_j}), (\ref{equation: decaying state, even N, x_j}).
Here $y_l/x_1$ is pure imaginary and $x_l/x_1$ is real. Hence, $\vec{\psi}^\alpha$ is not an eigenstate of $\mathcal{P}$ since the prefactor to $y_l$ is finite, i.e., $E_\pm \neq 0$. We conclude that in a finite Kitaev chain with even number of sites Majorana zero modes emerge only at the Kitaev points for $\mu=0$, since the states in Eqs. (\ref{equation: Delta = t, N even})-(\ref{equation: Delta = -t, N even}) are eigenstates of $\mathcal{P}$ as well. In the situation of odd $N$ and $\mu=0$, the eigenstates given by their elements in Eqs. (\ref{equation: zero mode on alpha N odd}), (\ref{equation: zero mode on beta N odd}) 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 $\xi$, while the case of full localisation is only reached at the the Kitaev points, where the MZM turn into the states given by Eqs. (\ref{equation: zero mode on alpha, odd N, Kitaev point plus})-(\ref{equation: zero mode on beta, odd N, Kitaev point minus}).
\section{Results for the spectrum and eigenstates at finite $\mu$}
\label{section: influence of non zero chemical potential}
\subsection{Spectrum}
The last missing situation is to consider a finite chemical potential $\mu$.
For this purpose we use  the so called chiral basis $\hat{\Psi}_c\mathrel{\mathop:}=\left(\textcolor{blue}{\gamma_1^A},\,\textcolor{blue}{\gamma_2^A},\,\ldots,\,\textcolor{blue}{\gamma_N^A},\,\textcolor{orange}{\gamma_1^B},\,\textcolor{orange}{\gamma_2^B},\,\ldots,\,\textcolor{orange}{\gamma_N^B}\right)^\mathrm{T}$. The Kitaev Hamiltonian transforms via $\hat{H}_\mathrm{KC}\,=\,\frac{1}{2}\,\hat{\Psi}^\dagger_c\,\mathcal{H}_c\,\hat{\Psi}_c,$
into a block off-diagonal matrix
\begin{align}\label{equation: Kitaev Hamiltonian in chiral basis}
	\mathcal{H}_c\,=\,\left[\begin{matrix}
	0_{N\times N} & h\\
	h^\dagger & 0_{N\times N}
\end{matrix}			\right],
\end{align}
because there are no $\gamma_j^A\,\gamma_i^A$ ($\gamma_j^B\,\gamma_i^B$) contributions in Eq. (\ref{equation: Kitaev Hamilt. in SSH basis, even N}). The $N \times N$ matrix $h$ is tridiagonal
\begin{align}\label{equation: matrix h in chiral basis}
	h=\left[\begin{matrix}
		-i\mu & a & \\
		-b & -i \mu & a\\
		& -b & -i \mu & a\\
		& & \ddots & \ddots & \ddots \\
		&&& -b & -i \mu & a\\
		&&&& -b & -i \mu
	\end{matrix}\right],
\end{align}
since the Kitaev Hamiltonian contains only nearest neighbour hoppings. Then the characteristic polynomial is \cite{Silvester}
\begin{align}\label{equation: characteristic polynomial via chiral basis}
	P_\lambda\left(\mathcal{H}_c\right)\,=\,\text{det}\left(\lambda^2\,\mathbbm{1}_N\,-\,h\,h^\dagger\right),
\end{align}
where, however, $h$ and $h^\dagger$ do not commute except for $t=0$ or $\Delta=0$. Thus, such matrices cannot be diagonalised simultaneously. Nevertheless the eigenvalues $\eta_j$ ($\eta_j^*$) of $h$ ($h^\dagger$) are easily derived e.g. following Ref.~[\onlinecite{Kouachi}]. We find
\begin{align}\label{equation: eigenvalues of h,h^dagger}
	\eta_j\,=\,-i\mu\,+\,2\,\sqrt{\Delta^2-t^2}\,\cos\left(\frac{j\pi}{N+1}\right),\quad j\,=\,1,\,\ldots,\,N,
\end{align}
independent of whether $\Delta\ge t$ or $t>\Delta$.
\subsubsection{Condition for zero energy modes}
Eq. (\ref{equation: eigenvalues of h,h^dagger}) immediately yields the criterion for hosting zero energy modes. According to Eq. (\ref{equation: Kitaev Hamiltonian in chiral basis}), we  have
	%
	%
	\begin{align}\label{equation: determinant of the Kitaev chain in terms of h and h*}
	\mathrm{det}(\mathcal{H}_c)\,=\,\mathrm{det}(h)\,\mathrm{det}(h^\dagger)\,=\, \vert \mathrm{det}(h)\vert^2,
	\end{align}
	%
	%
	and we need only to focus on $\mathrm{det}(h)$.
	%

	%
	If a single eigenvalue $\eta_j$ of $h$ is zero then  $\mathrm{det}(h)$ vanishes. Thus, for a zero energy mode the chemical potential must satisfy
	%
	%
	\begin{align}\label{equation: eigenvalue zero condition kouachi}
	\mu_n\,\,=\,2\,\sqrt{t^2-\Delta^2}\,\cos\left(\frac{n\pi}{N+1}\right),\quad n = 1,\,\ldots,\,N .
	\end{align}
	%
	%
	%
	%
 Obviously, Eq. (\ref{equation: eigenvalue zero condition kouachi}) cannot be satisfied for generic values of $t^2 -\Delta^2$, because all other quantities are real. The only possibility is $t^2\ge \Delta^2$. There is only one exception for odd $N$, because the value $n = (N+1)/2$ leads to $\mu = 0$ in Eq. (\ref{equation: eigenvalue zero condition kouachi}) for all values of $t$ and $\Delta$, in agreement with our results of section \ref{section: spectral analysis}. This result is exact and confirms findings from Ref.~[\onlinecite{Kao-2014}, \onlinecite{Hegde-2015}]; further it improves a similar but approximate condition on the chemical potential discussed by Zvyagin in Ref.~[\onlinecite{Zvyagin-2015}].
	\begin{figure}[th]
	\centering

		\includegraphics[width=\columnwidth]{figure9.pdf}
		\caption{Existence of zero energy solutions for growing system sizes. The red lines mark the boundaries of the topological phase diagram. Zero energy solutions correspond to the blue curves, along which the determinant of the Kitaev Hamiltonian vanishes as a function of $\mu$, $t$ and $\Delta$. a) The shown situation is for a small even $N$, yielding a zero determinant on the horizontal axis only at the Kitaev points $t/\Delta = \pm 1$. Each blue curve departs from one of these two points. b) The situation of small odd $N$ is similar to the even one, but the entire $\mu=0$-axis is now included. c) The solutions $\mu_n$ become dense for larger even $N$, but one sees already the filling of the non trivial phase for $N \rightarrow \infty$. d) Large even $N$ behave similar to large odd $N$, but the latter still include the entire horizontal axis.
	}
		\label{figure: finite size effects and the topological phase diagramm}
	\end{figure}
	%

	An illustration of these discrete solutions $\mu_n$, which we dub "Majorana lines", is shown in Fig. \ref{figure: finite size effects and the topological phase diagramm}. All paths contain the Kitaev points at $\mu = 0$ and $t = \pm \Delta$. 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 $\pi$.

	For growing number of sites $N$, the density of solutions increases. In the limit $N\rightarrow \infty$, $\theta_n= n\pi/N+1$ takes all values in $[0,\,\pi]$ and the entire area between $\mu  = \pm 2\,\sqrt{t^2-\Delta^2}$ for $t^2\ge \Delta^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.

\subsubsection{The complete spectrum of the finite Kitaev chain}

To proceed we transform Eq. (\ref{equation: characteristic polynomial via chiral basis}) into an eigenvector problem for $h h^\dagger$,
\begin{align}\label{equation: eigenvector problem of hhdagger}
	hh^\dagger \vec{v} \,=\,\lambda^2 \vec{v},
\end{align}
where we defined $\vec{v} = \left(\xi_1,\,\xi_2\,\ldots,\,\xi_N\right)^\mathrm{T}$.
Notice that we are not really interested in the eigenvector $\vec{v}$ here; we simply use its entries as dummy variables to release a structure hidden in the product of $h$ and $h^\dagger$. The elements of $h$,
\begin{align*}
	h_{n,m}\,=\, -i\mu\,\delta_{n,m}\,+\,a\, \delta_{n,m+1}\,-\,b\,\delta_{n+1,m},
\end{align*}
where $n,\,m \,= \,1,\ldots,N$, allow us to calculate the product $hh^\dagger$ entry wise
\begin{align*}
	\left(hh^\dagger\right)_{n,m}\,&=\,\delta_{n,m}\,\left[\mu^2\,-\,a^2\,\left(1-\delta_{n,N}\right)\,-\,b^2\,\left(1-\delta_{n,1}\right)
		\right]\\
		&\quad +\,i\mu\, \left(a-b\right)\,\left[\,\delta_{n,m+1}\,+\, \delta_{n+1,m}\right]\\
		&\quad +\,ab \,\left(\delta_{n,m+2}\,+ \,\delta_{n+2,m}\right).
\end{align*}
\begin{figure}[t]
	\centering

		\includegraphics[width = 1\columnwidth]{figure8.pdf}
	\caption{Spectrum of the Kitaev molecule with four sites and $\mu\neq0$. The green line follows the excitation spectrum from Eq. (\ref{infinite chain: excitation spectrum}) and the horizontal lines are the numerical eigenvalues of the Kitaev chain. The momenta $k_{3,\,2,\,1}$ are the proper wavevectors for $\mu \neq 0$ calculated from the full quantization rule, see Eqs.  (\ref{equation: relation k_1,2 via sum of cos}), (\ref{equation: Kitaev chain full quantization rule}), (\ref{equation: definition of F for the quantization rule}). The dashed, light grey lines represent the wavevectors taken from the $\mu=0$ case to highlight the difference. The chemical potential $\mu$ is obviously changing the quantization of a finite chain. The chosen parameters are $t=4\,$meV, $\Delta=1.5\,$meV and $\mu=3\,$meV and lead to the numerical energies $\left[ 0.43,\,  4.034,\,  6.068,\,  9.603\right]$ (in meV). The value of $k_{3,\,2,\,1} d$ is approximately $\left[0.6360,\, 1.2753,\,1.6086\right]$. }
	\label{figure: mu influences the quantization}
\end{figure}
Thus, importantly, Eq. (\ref{equation: eigenvector problem of hhdagger}) reveals a recursion formula
\begin{align}\label{equation: tetranacci recursion formula}
	\xi_{j+2}\,&=\,\frac{\lambda^2+a^2+b^2-\mu^2}{ab}\,\xi_j-\xi_{j-2}\notag\\
		&\quad -i\mu\left(\frac{a-b}{ab}\right)\,\left(\xi_{j+1} + \xi_{j-1}\right),
\end{align}
for the components of $\vec{v}$. The entries $\xi$ are a generalisation of the Fibonacci polynomials $\zeta_j$ from Eq. (\ref{equation: recursion formula fibonacci polynomials}), to which they reduce for $\mu= 0$, and may be called Tetranacci polynomials\cite{McLaughlin-1979, Waddill-1992}.
Further,  we find the open boundary conditions from Eq. (\ref{equation: eigenvector problem of hhdagger}) to be
\begin{align}\label{equation: open boundary condition in general}
	\xi_0 \, =\, \xi_{N+1}\,=\,b\,\xi_{N+2}\,-\,a\, \xi_N\,=\,b\,\xi_1\,-\,a \,\xi_{-1}\,=\,0,
\end{align}
where we used Eq. (\ref{equation: tetranacci recursion formula}) for simplifications.

Appendix \ref{appendix: spectrum finite u} contains the description of how to deal with those polynomials, the boundary conditions and further the connection of Eq. (\ref{equation: tetranacci recursion formula}) to Kitaev's bulk spectrum $\lambda = E_\pm(k)$ in Eq. (\ref{infinite chain: excitation spectrum}). Essentially one has to use similar techniques as it was  done for the Fibonacci polynomials, where now the power law ansatz $\xi_j\propto r^j$ leads to a characteristic equation for $r$ of order four. Thus, we find in total four linearly independent fundamental solutions $r_{\pm 1\,, \pm 2}$, which can be expressed in terms of two complex wavevectors denoted by $\kappa_{1,2}$ through the equality
\begin{align}
	r_{\pm j}\,=\,e^{\pm i \kappa_j },\quad j= 1,2.
\end{align}
These wavevectors are not independent, but coupled via
\begin{align}\label{equation: relation k_1,2 via sum of cos}
	\cos(\kappa_1)+ \cos(\kappa_2)\,=\,-\frac{\mu t}{t^2-\Delta^2},\quad \forall \,t,\,\Delta,\, \mu\in\mathbbm{R}.
\end{align}
\begin{widetext}

\begin{center}
	\centering
	\begin{table}
	\caption{Overview of the quantization rule for the wave vectors of the finite Kitaev chain in different scenarios. The wavevectors $k$, $k_n$, $\kappa_{1,2}$ used together with Eq. (\ref{infinite chain: excitation spectrum}) and $q$ with Eq. \eqref{equation: spectrum: mu=0, even N, edges} yield the correct finite system energies and $k,\,k_n,\,q\in\mathbbm{R}$, $\kappa_{1,2}\in\mathbbm{C}$. Notice: $n = 1,\,\ldots,\,N$ and $m = \frac{t}{\Delta}$ ($m = \frac{\Delta}{t}$) for $\vert t \vert > \vert \Delta\vert$ ($\vert \Delta \vert > \vert t\vert$).}
	\begin{ruledtabular}

	\begin{longtable}{p{2.4 cm} p{7 cm} p{4.2 cm} p{2.2 cm} p{1.4 cm} }

	\hline\hline \\

	Requirements & Quantisation rule &  Zero modes & Equation for \newline eigenstate \newline elements  & Majorana\newline character\\
	\\
	\hline\\
	$\Delta = 0:$		 & $k_n d = \frac{n\pi}{N+1}$						& Yes, if for some $n$: \newline $\mu=\mu_n = 2t \cos\left(\frac{n\pi}{N+1}\right) $ &		& No \\
	\\
	$t = 0:$ 			 & $k_n d = \frac{\pi}{2}+\frac{n\pi}{N+1}$						& No & 						&	No	\\
	\\ \\ \\
	$\mu = 0,$ $N$ odd:  & $k_n d = \frac{n\pi}{N+1},\quad n\neq \frac{N+1}{2}$ \newline\newline
						   $qd =\mathrm{arctanh}\left(1/\vert m \vert\right)$	 		& No \newline\newline Yes & (\ref{equation: N odd, entries of the eigenstates})\,-\,(\ref{equation: function T odd })\newline\newline
						   (\ref{equation: zero mode on alpha N odd}), (\ref{equation: zero mode on beta N odd})\newline (\ref{equation: zero mode on alpha, odd N, Kitaev point plus})\,-\,(\ref{equation: zero mode on beta, odd N, Kitaev point minus})	&No	\newline\newline Yes\newline Yes	\\
	\\ \\ \\
	$\mu = 0,$ $N$ even: & $\tan\left[kd \left(N+1\right)\right]\,=\,\mp\frac{\Delta}{t}\tan\left(kd\right)$\newline\newline
						  $\tanh\left[qd \left(N+1\right)\right]\,=\,\mp m\tanh\left(qd\right)$ & No \newline\newline only if $\Delta = \pm t$\newline otherwise	& (\ref{equation: decaying state, even N, y_j}), (\ref{equation: decaying state, even N, x_j})\newline\newline (\ref{equation: Delta = t, N even})\,-\,(\ref{equation: Delta = -t, N even})\newline (\ref{equation: decaying state, even N, y_j}), (\ref{equation: decaying state, even N, x_j}) & No	\newline\newline Yes \newline No	\\
	\\ \\ \\
	$t,\,\Delta,\,\mu\in\mathbbm{R}$ 				&  $\displaystyle\frac{\sin^2\left[\frac{\kappa_1+\kappa_2}{2}  \left(N+1\right)\right]}{\sin^2\left[\frac{\kappa_1-\kappa_2}{2}  \left(N+1\right)\right]} = \frac{1+\left(\frac{\Delta}{t}\right)^2\,\cot^2\left(\frac{\kappa_1-\kappa_2}{2}\right) }{1+\left(\frac{\Delta}{t}\right)^2\,\cot^2\left(\frac{\kappa_1+\kappa_2}{2}\right)}$ &  only for $\mu=\mu_n\in \mathbbm{R}$ and \newline \makebox{$\mu_n = 2\sqrt{t^2-\Delta^2} \cos\left(\frac{n\pi}{N+1}\right) $}	&\eqref{equation: closed formula for all Tetranacci xi} \eqref{equation: entries of u}\newline(\ref{equation: MZM x_l}), (\ref{equation: MZM ycal_l}) \newline (\ref{equation: MZM xcal_l}), (\ref{equation: MZM y_l})  & No \newline Yes\newline Yes\\
	\\
	\hline\hline
	\end{longtable}

	\end{ruledtabular}
	\label{table: all quantisiation (sub) rules}
	\end{table}
\end{center}

\end{widetext}

For $\mu=0$ we can recover from Eq. (\ref{equation: relation k_1,2 via sum of cos}) our previous results, whereby one has only pure real $(k)$ or pure imaginary $(iq)$ wavevectors\footnote{In fact $\mu=0$ supports complex wavevectors too, but their real part has to be zero or $\pi/2$, i.e. one has to use $iq$ or $\frac{\pi}{2}+iq$.}.
Further, Eqs. (\ref{infinite chain: excitation spectrum}) and (\ref{equation: relation k_1,2 via sum of cos}) yield
\begin{align*}
	E_\pm (\kappa_1) = E_\pm(\kappa_2).
\end{align*}
The linearity of the recursion formula Eq. (\ref{equation: tetranacci recursion formula}) states that the superposition of all four fundamental solutions is the general form of $\xi_j$. Since the boundary conditions translate into a homogeneous system of four coupled equations and a trivial solution for $\xi_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
\begin{align}\label{equation: Kitaev chain full quantization rule}
	F(\kappa_1, \kappa_2) = F(\kappa_1, -\kappa_2),
\end{align}
where we introduced the function $F(\kappa_1, \kappa_2)$ as
\begin{align}\label{equation: definition of F for the quantization rule}
	F(\kappa_1, \kappa_2)&\,=\,\sin^2\left[\frac{\kappa_1+\kappa_2}{2}\,\left(N+1\right)\right]\times\notag\\
	&\quad \times \left[1+\left(\frac{\Delta}{t}\right)^2\,\cot^2\left(\frac{\kappa_1+\kappa_2}{2}\right)\right].
\end{align}

Similar quantization conditions are known for an open X-Y spin chain in transverse field\cite{Loginov-1997}.
Notice that the quantization rule is symmetric with respect to $\kappa_{1,2}$. \mbox{Table \ref{table: all quantisiation (sub) rules}} gives an overview of the quantization rules for different parameter settings $\left(\Delta,\,t,\,\mu\right)$. The bulk eigenvalues of a finite Kitaev chain with four sites and $\mu\neq 0 $ are shown in Fig. \ref{figure: mu influences the quantization}.

The previous relations open another route to finding the condition   leading to modes with exact zero energy. A convenient form of Eq. (\ref{equation: relation k_1,2 via sum of cos}) is
\begin{align}\label{equation: connection of kappa1pmkappa2}
	\cos\left(\frac{\kappa_1+\kappa_2}{2}\right)\, \cos\left(\frac{\kappa_1-\kappa_2}{2}\right)\,=\,-\frac{1}{2}\,\frac{\mu t}{t^2-\Delta^2},
\end{align}
and the dispersion relation can be transformed into
\begin{align}\label{equation: dispersion relation with kappas}
	E^2&=\frac{1}{\cos^2\left(\frac{\kappa_1\pm\kappa_2}{2}\right)}
	\left[4(t^2-\Delta^2)\cos^2\left(\frac{\kappa_1\pm\kappa_2}{2}\right)-\mu^2  \right]\times \notag	\\
	&\quad \times\left[\frac{t^2}{t^2-\Delta^2}-\cos^2\left(\frac{\kappa_1\pm\kappa_2}{2}\right)\right].
\end{align}
Both combinations $\kappa_1\pm\kappa_2$ yield the same energy, due to Eq. (\ref{equation: connection of kappa1pmkappa2}). If one of the brackets in Eq. (\ref{equation: dispersion relation with kappas}) vanishes for $\kappa_1+\kappa_2$ ($\kappa_1-\kappa_2$), the second one does so for $\kappa_1-\kappa_2$ ($\kappa_1+\kappa_2$ ) too. Hence, zero energy is achieved exactly if
\begin{align}\label{equation: cot part in the quantization rule is zero}
	1+\left(\frac{\Delta}{t}\right)^2\,\cot^2\left(\frac{\kappa_1\pm\kappa_2}{2}\right)\,=\,0.
\end{align}
This  puts restrictions on $\Delta,\,t,\,\mu$. Together with  the quantization rule in Eq. (\ref{equation: Kitaev chain full quantization rule}), this ultimately leads to (\ref{equation: eigenvalue zero condition kouachi}) and  to the condition $2\sqrt{t^2-\Delta^2} < \vert \mu \vert < 2 \vert t \vert$, 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\rightarrow \infty$ the difference between even and odd $N$ vanishes and the part of the $\mu=0$ axis between $t=-\Delta$ and $t=\Delta$ for even $N$ leads to zero energy states too, in virtue of Eq. (\ref{equation: spectrum: mu=0, even N, edges, cosh/cosh version}). Analogously, the area around the origin in Fig. (\ref{figure: finite size effects and the topological phase diagramm}), defined by $2\sqrt{t^2-\Delta^2}>\vert \mu \vert $ and  $\mu < 2 \vert t \vert$ with $\mu \neq 0$ does not  support zero energy modes for all finite $N$.

 Instead, this area contains solutions with exponentially small energies, see Eq. (\ref{equation: dispersion relation with kappas}), which become zero exclusively in the limit $N\rightarrow \infty$. The wavevectors obey $\kappa_1 = i q_1$, $\kappa_2 = \pi + i q_2$ with real $q_{1,2}$ for $\Delta^2>t^2$, and  $\kappa_1 = i q_1$, $\kappa_2 =  i q_2$ otherwise, which follows from Eq. (\ref{equation: Kitaev chain full quantization rule}) after some manipulations; see also Ref.~[\onlinecite{Loginov-1997}]. Thus, the entire non trivial phase hosts zero energy solutions for $N\rightarrow\infty$.
\subsection{Eigenstates}
The calculation of an arbitrary wave function of the Kitaev Hamiltonian without any restriction on $t$, $\Delta$, $\mu$ is performed at best in the chiral basis yielding the block off-diagonal structure in Eq. \eqref{equation: Kitaev Hamiltonian in chiral basis}. A suitable starting point is to consider a vector $\vec{w}$ in the following form
\begin{align*}
	\vec{w}\,=\,\left(
	\begin{matrix}
	\vec{v}\\
	\vec{u}
	\end{matrix}
	\right)
\end{align*}
with $\vec{v}=\left(\xi_1\,\ldots,\,\xi_N\right)$, $\vec{u}=\left(\sigma_1\,\ldots,\,\sigma_N\right)$. Solving for an eigenstates with eigenvalue $\lambda$ demands on $\vec{v},\,\vec{u}$
\begin{align}
    \label{eq:hu-is-lambdav}
	h\,\vec{u}\,&=\,\lambda\,\vec{v},\\
	\label{equation: u follows from v}
	h^\dagger\,\vec{v}\,&=\,\lambda\,\vec{u},
\end{align}
with $h$ from Eq. \eqref{equation: matrix h in chiral basis}. Thus,
$	h\,h^\dagger \vec{v}\,=\,\lambda^2\,\vec{v}$,
and we recover Eq. \eqref{equation: eigenvector problem of hhdagger} and the entries of $\vec{v}$ obey Eq. \eqref{equation: tetranacci recursion formula} again. In appendix \ref{appendix: section: The closed formula of Tetranacci polynomials} we derived the closed formula for $\xi_j$, namely
\begin{align}\label{equation: closed formula for all Tetranacci xi}
	\xi_j\,=\,\sum\limits_{i=-2}^1\,\xi_i \,X_i(j), \quad j \in \mathbbm{Z}
\end{align}
where $\xi_{-2},\,\ldots \xi_1$ are the initial values of the polynomial sequence dependent on the boundary conditions, and $X_i(j)$ inherit the selective property
\begin{align}
	X_i(j)\,=\,\delta_{i,j}\qquad \mathrm{for~only\,} i,\,j=-2,\,\ldots,\,1.
\end{align}
That these functions $X_i(j)$ exist and that they indeed satisfy Eq. \eqref{equation: closed formula for all Tetranacci xi} for arbitrary values of $j$ is discussed in Appendix \ref{appendix: section: The closed formula of Tetranacci polynomials}.

The remaining task is to obtain the initial values $\xi_{-2},\,\ldots \xi_1$,
 which follow from the open boundary conditions  Eq. \ref{equation: open boundary condition in general}. Further, one has one free degree of freedom, which we  to choose to be the  entry $\xi_1$ of $\vec{v}$. In total our initial values are $\xi_1,\xi_0=0,\xi_{-1}=b \,\xi_1/a$ and $\xi_{-2}$ follows from $\xi_{N+1}=0$ and Eq. \eqref{equation: closed formula for all Tetranacci xi},
\begin{align*}
	\xi_{-2}\,=\,-\xi_1\,\frac{a\,X_1(N+1)+b \,X_{-1}(N+1)}{a\,X_{-2}(N+1)}.
\end{align*}
Demanding further $b\,\xi_1\,-\,a \,\xi_{-1}\,=\,0$ quantizes  the momenta $\kappa_{1,2}$ and in turn $\lambda=E_\pm(\kappa_{1,2})$, according to Eq. \eqref{equation: Kitaev chain full quantization rule} (the relation between $X_j$ and $\kappa_{1,2}$ is discussed in the Appendix \ref{appendix: section: The closed formula of Tetranacci polynomials}). Notice that the form given by  Eqs. \eqref{equation: closed formula for all Tetranacci xi}, \eqref{equation: entries of u} (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 $\kappa_{1,2}$.

The second part of the eigenstate $\vec{w}$, i.e. $\vec{u}$,  follows in principle from Eq. \eqref{equation: u follows from v}. A simpler and faster way is to consider
 $	h^\dagger\,h\,\vec{u}\,=\,\lambda^2\,\vec{u}$.
A comparison of $h$ and $h^\dagger$ reveals that they transform into each other by exchanging $a$ and $b$ and switching $\mu$ into $-\mu$. Consequently, the structure of the entries of $\vec{u}$ follow essentially from the ones of $\vec{v}$. Thus, we find that the $\sigma_i$ obey
Eq. \eqref{equation: tetranacci recursion formula} too, yielding
\begin{align}\label{equation: entries of u}
	\sigma_j\,=\,\sum\limits_{i=-2}^1\,\sigma_i \,X_i(j), \quad j \in \mathbbm{Z}
\end{align}
with the same functions $X_i(j)$. The boundary condition on $\vec{u}$ reads
\begin{align*}
	\sigma_0 \, =\, \sigma_{N+1}\,=\,a\,\sigma_{N+2}\,-\,b\, \sigma_N\,=\,a\,\sigma_1\,-\,b \,\sigma_{-1}\,=\,0.
\end{align*}
Proceeding as we did for $\vec{v}$ yields $\sigma_0\,=\,0$, $\sigma_{-1}=a \sigma_1/b$ and
\begin{align}
	\sigma_{-2}\,=\,-\sigma_{1}\,\frac{b\,X_1 (N+1)-\,a\,X_{-1}(N+1)}{b\,X_{-2}(N+1)},
\end{align}
where $\sigma_1$ is fixed by the first line of Eq. \eqref{equation: u follows from v}
\begin{align}
\label{eq:sigma1-with-lambda}
	\sigma_1\,=\,\frac{i\mu\,\xi_1\,+\,b\,\xi_2}{\lambda}.
\end{align}
The last open boundary condition $a\,\sigma_{N+2}\,-\,b \,\sigma_{N}\,=\,0$ is satisfied for $\lambda=E_\pm(\kappa_{1,2})$ and thus already by the construction of $\vec{v}$. We notice that we assumed $\lambda\neq 0$ in order to obtain $\sigma_1$; the  case $\lambda=0$ is discussed in Sec. \ref{section: MZM paths in the topological phase diagram for finite N} below.

In the limit of $\mu=0$ on $\vec{w}$ it holds that
\begin{align*}
	X_{-2}(2l+1)\vert_{\mu=0} \,&=\,0,\\
	X_{0}(2l+1)\vert_{\mu=0} \,&=\,0,\\
	X_{-1}(2l)\vert_{\mu=0} \,&=\,0,\\
	X_{1}(2l)\vert_{\mu=0} \,&=\,0,
\end{align*}
for all values of $l$; thus only two initial values for $\xi_j$ and two for $\sigma_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 $\xi_j$ ($\sigma_j$) are real (pure imaginary) if $\xi_1$ is chosen to be real. Thus, the corresponding operators $\psi_{\vec{w}}$  $(\psi^\dagger_{\vec{w}})$ can never square to $1/2$.

\section{MZM eigenvectors at finite $\mu$}
\label{section: MZM paths in the topological phase diagram for finite N}

The technique outlined above (in particular, Eq.~\eqref{eq:sigma1-with-lambda}) cannot be used directly for exact zero energy modes, because then $\lambda = 0$ in Eqs.~\eqref{eq:hu-is-lambdav},~\eqref{equation: u follows from v}.
In this section we demonstrate the Majorana nature of the zero energy solutions satisfying Eq. (\ref{equation: eigenvalue zero condition kouachi}), and we give the explicit form of the associated MZM using a different technique, which (similar to the Chebyshev polynomials method in Ref.~[\onlinecite{Kawabata-2017}]) requires only the use of Fibonacci, not Tetranacci polynomials. This simplification is caused by the fact that setting $\lambda=0$ decouples the two Majorana sublattices, while setting $\mu=0$ decouples the two SSH-like chains.
\begin{figure}[htb]\centering

	\includegraphics[width=\columnwidth]{figure10.pdf}
	\caption{Illustration of the sawtooth pattern of $\vec{\psi}_A$. The real space position of the entries $x_l$ (${\scriptstyle{\mathcal{Y}}}_l$) of $\vec{\psi}_A$ is
	at $j=2l-1$ ($j=2l$) on chain $\alpha $ ($\beta$) and marked by the blue (light blue) spheres. The blue line connects these entries as guide to the eye.
		}
	\label{Figure: sawtooth pattern of  psi_A}
\end{figure}

We use the SSH-like description Eq. (\ref{equation: Kitaev Hamiltonian/ matrix in SSH basis}) of the Kitaev chain where $\mu\neq 0 $ couples both chains together. Consequently an eigenstate $\vec{\psi}\,=\,\left(\vec{v}_\alpha,\,\vec{v}_\beta\right)^\mathrm{T}$ has in general no zero entries and we use the same notation for the components of $\vec{v}_\alpha$, $\vec{v}_\beta$ as in the sections \ref{subsection: N even eigenvectors} and \ref{subsection: N odd eigenvectors}.

The zero energy values are twice degenerated, as one can see from Eq. (\ref{equation: determinant of the Kitaev chain in terms of h and h*}), and the associated zero modes are connected by the chiral symmetry $\mathcal{C}$. Thus, we get zero energy states by superposition $\vec{\psi}_{A,B} \,\mathrel{\mathop:}= \left(\vec{\psi}\pm \mathcal{C}\vec{\psi}\right)/2$ too. The chiral symmetry Eq. (\ref{equation: chiral symmetry even N}), contains an alternating pattern of $\pm 1$, such that $\vec{\psi}_A$ $(\vec{\psi}_B)$ includes only non zero entries on the Majorana sublattice $A$ ($B$). Hence, $\vec{\psi}_A$
($\vec{\psi}_B$) contains only $x_l$ (${\scriptstyle{\mathcal{X}}}_l$) and ${\scriptstyle{\mathcal{Y}}}_j$ ($y_j$) terms and the last component depends on whether $N$ is odd or even. In the latter case we have
\begin{align*}
	\vec{\psi}_A\,&=\,\left(x_1,\,0,\,x_2,\,0,\ldots,x_{\frac{N}{2}},0\,\left\vert\,0,{\scriptstyle{\mathcal{Y}}}_1,\,0,{\scriptstyle{\mathcal{Y}}}_2,\ldots,0,\,{\scriptstyle{\mathcal{Y}}}_{\frac{N}{2}}\right.\right)^\mathrm{T},\\
	\vec{\psi}_B\,&=\,\left(0,\,y_1,\,0,\,y_2,\ldots,0,y_{\frac{N}{2}}\left\vert\,{\scriptstyle{\mathcal{X}}}_1,\,0,{\scriptstyle{\mathcal{X}}}_2,\,0\ldots,\,{\scriptstyle{\mathcal{X}}}_{\frac{N}{2}},\,0\right.\right)^\mathrm{T}.
\end{align*}
The form of the odd $N$ eigenvectors is quite similar, see Eqs. (\ref{equation: appendix: vec_+ for odd N}), (\ref{equation: appendix: vec_- for odd N}).

The composition of $\vec{\psi}_A$ is illustrated in Fig. \ref{Figure: sawtooth pattern of  psi_A}, where its entries are shown to form a sawtooth like pattern, following the action of $\gamma_j^A$ on both SSH-like chains.

The full calculation is given in appendix \ref{appendix: zero energy states}. We focus here on $\vec{\psi}_A$ exclusively, because the $\vec{\psi}_B$ components follow essentially from $\vec{\psi}_A$ by exchanging $a$ and $b$ and replacing $i\mu$ by $-i\mu$. The chemical potential has still to obey Eq. (\ref{equation: eigenvalue zero condition kouachi}).

The components of the zero mode $\vec{\psi}_A$ have to satisfy $(l = 1,\,\ldots,\,N/2-1)$
\begin{align}\label{equation: recursion formula x_l}
	b\,x_{l+1}\,-\,a\,x_l\,+\,i\mu\,{\scriptstyle{\mathcal{Y}}}_l\,&=\,0,\\
	b\,{\scriptstyle{\mathcal{Y}}}_{l+1}\,-\,a\,{\scriptstyle{\mathcal{Y}}}_l\,+i\mu\,x_{l+1}\,&=\,0,
\end{align}
for even $N$, and the open boundary conditions are
\begin{align*}
	{\scriptstyle{\mathcal{Y}}}_0\,=\,y_0\,=\,{\scriptstyle{\mathcal{X}}}_{\frac{N}{2}+1} \,=\,x_{\frac{N}{2}+1} \,=\,0.
\end{align*}
The situation for the entries of $\vec{\psi}_A$ for odd $N$ is similar
\begin{align}
	b\,x_{j+1}\,-\,a\,x_j\,+\,i\mu\,{\scriptstyle{\mathcal{Y}}}_j\,=\,0,\\
	\label{equation: recursion formula ycal_l}
	b\,{\scriptstyle{\mathcal{Y}}}_{i+1}\,-\,a\,{\scriptstyle{\mathcal{Y}}}_i\,+i\mu\,x_{i+1}\,=\,0,
\end{align}
where $j = 1,\ldots,\,(N-1)/2$, $i = 1,\ldots,\,(N-3)/2$. The open boundary condition changes to
\begin{align*}
	{\scriptstyle{\mathcal{Y}}}_0\,=\,y_0\,=\,{\scriptstyle{\mathcal{Y}}}_{\frac{N+1}{2}} \,=\,y_{\frac{N+1}{2}} \,=\,0.
\end{align*}
Solving these recursive formulas leads in both cases to
\begin{align}\label{equation: MZM x_l}
	x_l\,&=\,x_1\,\frac{\sin[\theta_n\,(2l-1)]}{\sin(\theta_n)}\,\left(-\frac{a}{b}\right)^{l-1},
\end{align}
with $\theta_n = n\pi/(N+1)$, $n =1,\,\ldots,\,N$ and
\begin{align}\label{equation: MZM ycal_l}
	{\scriptstyle{\mathcal{Y}}}_l\,&=\,-x_1\,\mathrm{sgn}(t+\Delta)\frac{\sin(2\,\theta_n\,l)}{\sin(\theta_n)}\,\left(-\frac{a}{b}\right)^{\frac{2l-1}{2}},
\end{align}
where $x_1$ is a free parameter and $-a/b\ge 0$ due to $t^2\ge \Delta^2$. Recalling that $a = i \left(\Delta -t\right)$,  $b = i \left(\Delta +t\right)$, Eqs. (\ref{equation: MZM x_l}), (\ref{equation: MZM ycal_l}) predict an oscillatory exponential decay of the
coefficients $x_l$, ${\scriptstyle{\mathcal{Y}}}_l$. For example
\begin{align}
	x_l\,&=\,x_1\,\frac{\sin[\theta_n\,(2l-1)]}{\sin(\theta_n)}\,e^{-\left(l-1\right) d/\xi},
\end{align}
where the decay length is defined by $\xi= d/ \left\vert\ln \left(\frac{t-\Delta}{t+\Delta}\right)\right\vert$, for $t>\Delta>0$. Summarizing: the zero energy modes $\vec{\psi}_{A,B}$ look like small or strong suppressed standing waves with $n-1$ nodes for $n=1,...,N_\mathrm{max}$ and $N_\mathrm{max}  = N/2$  ($N_\mathrm{max} = N-1/2$) for even (odd) $N$. The expressions for ${\scriptstyle{\mathcal{X}}}_l$ and $y_l$ are obtained in a similar way
\begin{align}\label{equation: MZM xcal_l}
	{\scriptstyle{\mathcal{X}}}_l\,&=\,{\scriptstyle{\mathcal{X}}}_1\,\frac{\sin[\theta_n\,(2l-1)]}{\sin(\theta_n)}\,\left(-\frac{b}{a}\right)^{l-1},\\
	\label{equation: MZM y_l}
	y_l\,&=\,-{\scriptstyle{\mathcal{X}}}_1\,\mathrm{sgn}(t-\Delta)\frac{\sin(2\,\theta_n\,l)}{\sin(\theta_n)}\,\left(-\frac{b}{a}\right)^{\frac{2l-1}{2}},
\end{align}
and ${\scriptstyle{\mathcal{X}}}_1$ can be freely chosen. The open boundary conditions for $l=0$ are satisfied by construction of ${\scriptstyle{\mathcal{Y}}}_l$ ($y_l$), while the remaining ones follow due to the structure of $\theta_n$.
\begin{figure}[ht]
	\centering

	\includegraphics[width=\columnwidth]{figure11.pdf}
	\caption{Majorana zero mode $\vec{\psi}_A$ for various parameter sets. The considered parameters are denoted by different symbols on the topological phase diagram. The dark (light) blue spheres follow $x_l/x_1$ (${\scriptstyle{\mathcal{Y}}}_l/x_1$) at position $j=2l-1$ ($j=2l$) from Eq. (\ref{equation: MZM x_l}) (\ref{equation: MZM ycal_l}). The decay length of the MZM increases for larger ratios of $t/\Delta$ for a fixed value of $\theta_n$, until the state is delocalised over the entire system. Lowering the chemical potential, e.g. following the vertical orange line, but keeping $t/\Delta$ fixed, changes the shape of the MZM's. Small decaying length and large enough chemical potentials leads to Majorana modes which have highest weight in the center of the chain.
}
\label{figure: Majorana zero modes}
\end{figure}

The zero mode $\vec{\psi}_A$ is shown in Fig. \ref{figure: Majorana zero modes} for a various range of parameters. For not too large ratios $ t/\Delta  >1$, the zero mode $\vec{\psi}_A$ is mostly localised at one end of the Kitaev chain and decays away from it in an oscillatory way. The eigenstate $\vec{\psi}_B$ is concentrated on the opposite end. While the oscillation depends on the chemical potential $\mu_n$ associated to the zero mode, according to Eq. (\ref{equation: eigenvalue zero condition kouachi}), the decay length is only set by the parameters $\Delta$ and $t$. Thus, as the ratio of $t/\Delta$ is increased, the zero energy mode gets more and more delocalized.

The zero energy states $\vec{\psi}_{A,B}$ are MZM's, since they are eigenstates of the particle hole operator $\mathcal{P}$ Eq. (\ref{equation: particle hole symmetry in real space in ssh basis}) for real or pure imaginary values of $x_1$, ${\scriptstyle{\mathcal{X}}}_1$. Further, the states $\vec{\psi}=\vec{\psi}_A + \vec{\psi}_B$ and $\mathcal{C} \vec{\psi}=\vec{\psi}_A - \vec{\psi}_B$ are MZM's too. On the other hand a fermionic state is constructed with $\vec{\psi}_\pm = \vec{\psi}_A \pm i\, \vec{\psi}_B$, similar to what was found in Sec. \ref{section: Eigenvectors and symmetries} for the $\mu=0$ case, or at the Kitaev points in Eqs. \eqref{equation: Missing electron a} and \eqref{equation: Missing electron b} in Sec. \ref{section 2 model}.

There are three limiting situations we would like to discuss: $t\rightarrow \pm \infty$, $N\rightarrow \infty$, 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 $\xi$. Because $-a/b \rightarrow 1$ for $t\rightarrow \pm\infty$, this implies also that $\xi\rightarrow \infty$ in that limit. Hence oscillations are less suppressed for large values of $t$, as illustrated in Fig. \ref{figure: Majorana zero modes}. Already a ratio of $t/\Delta \approx 100$ is enough to avoid a visible decay for $N\approx 20$. This effect can be found as long as $N$ is finite, but one has to consider larger values of the ratio $t/\Delta$.

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$ ${\scriptstyle{\mathcal{Y}}}_l$, ${\scriptstyle{\mathcal{X}}}_l$ and $y_l$ leads to significantly large or small values. Thus, the state $\vec{\psi}_A$ ($\vec{\psi}_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$, ${\scriptstyle{\mathcal{Y}}}_l$ only change their sign. For odd $N$ and for $\theta_n = \pi/2$, i.e. $\mu=0$, one recovers the result in the Eqs. (\ref{equation: zero mode on alpha N odd}) -  (\ref{equation: zero mode on beta N odd}).

\section{Numerical results and impact of disorder}
\label{section: 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.
\subsection{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 Eq. \eqref{equation: Kit. Hamiltonian, fermionic operators, realspace} and varying parameters, are shown in Fig.~\ref{figure: numerical data groundstate energy}. The phase diagram in Fig.~\ref{figure: numerical data groundstate energy} (a) is the numerical equivalent of that shown in Fig.~\ref{figure: Majorana zero modes}, but for a smaller range of $t$ and $\mu$. 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, $\mu \lesssim 2t$, all the boundary states in a finite system have finite energy, as shown in Fig.~\ref{figure: numerical data groundstate energy} (b). Figure ~\ref{figure: numerical data groundstate energy} (c) displays $E_0$ for fixed $t/\Delta = 17$, as a function of $\mu$ and $N$. The number of near-zero energy solutions increases linearly with $N$, according to Eq.~\eqref{equation: eigenvalue zero condition kouachi}.
\begin{figure}[ht]
	\centering

	\includegraphics[width=\columnwidth]{figure12.pdf}
	\caption{Numerical results for the energy $E_0$ of the lowest positive energy state as a function of $t/\Delta$, $\mu/\Delta$ and system size $N$. (a) $E_0$ as a function of $t$ and $\mu$ for $N=20$. The red line marks the boundary of the bulk topological phase. The $N/2$ dark lines coincide with the Majorana lines given by Eq. (\ref{equation: eigenvalue zero condition kouachi}). (b) Zoom into the neighbourhood of the Kitaev point at $\mu=0$ and $t=\Delta$, showing the absence of zero energy solutions for $\Delta<t$ in the nominally non-trivial phase. (c) $E_0$ as a function of $\mu$ and the system size $N$ for $t/\Delta = 17$. The red line marks $\mu=2t$, the boundary of the bulk topological phase. As $N$ increases, the number of $\mu$ yielding zero energy solutions also increases according to Eq. (\ref{equation: eigenvalue zero condition kouachi}), and the maximum energies of the bound states decrease. (d) The values of $E_0$ for the same set of parameters as in (c), projected onto the $N$-$E_0$ plane. The ground state energy $E_0$ shows a $\mu$ dependent, oscillatory behavior in the system size $N$, where the maximum energies follow with very good accuracy an $a_\mu \exp(-Nd/\xi)$ rule for sufficiently large values of $N$, with $\xi$ from Eq. (\ref{equation: decay length}) and the numerical prefactor $a_\mu$.}
	\label{figure: numerical data groundstate energy}
\end{figure}
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 Fig.~\ref{figure: Majorana zero modes}). The decay length $\xi$ of the in-gap eigenstates, defined in Eq.~\eqref{equation: decay length}, is determined by the ratio $t/\Delta$ 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 {\it maximum} energy of the boundary states that decreases as $E_{0,max} \propto \exp(-Nd/\xi)$ as $N$ is increased, as illustrated in Fig.~\ref{figure: numerical data groundstate energy} (d), in agreement with Eq. \eqref{equation: spectrum: mu=0, even N, edges, cosh/cosh version} and Ref. [\onlinecite{kitaev:physusp2001}, \onlinecite{Zvyagin-2015}, \onlinecite{Zeng-2019}]. The {\it minimum} energy of zero can be reached for any chain length, provided that the chemical potential is appropriately tuned.

\subsection{Topological protection against Anderson disorder}
\label{subsection: disorder}
\begin{figure}[b]
	\begin{center}
		\includegraphics[width=\columnwidth]{Figure13.pdf}
	\end{center}
	\caption{\label{figure: disorder}
		The lowest ($E_0$) and second lowest ($E_1$) energy eigenvalue of a Kitaev chain with $N=20$ sites as a function of $\mu$ for various ratios of $t/\Delta$. Thin lines correspond to the eigenvalues of a clean chain, thick lines to those of a chain with random on-site disorder $\varepsilon_i \in [-4\Delta,4\Delta]$. The red line marks the boundary of the topological phase.
	}
\end{figure}
 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 $\mu$ for $N=20$ and three different values of $t/\Delta$. The disorder was modeled as an on-site energy term $\varepsilon_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 Fig.~\ref{figure: disorder}. In all plots $W=4\Delta$, i.e. twice larger than the $\pm2\Delta$ gap at $\mu=0$. Each curve is an average over 100 disorder configurations.
\begin{figure}[hbt]
	\begin{center}
		\includegraphics[width=\columnwidth]{Figure14.pdf}
	\end{center}
	\caption{\label{figure: disorder-noaverage}
		The lowest energy eigenvalue $E_0$ of a Kitaev chain with $N=20$ sites as a function of $\mu$ for $t/\Delta = 17$. Thin lines correspond to the eigenvalues of a clean chain, thick lines to those of a chain with random on-site disorder $\varepsilon_i \in [-W,W]$, for two disorder realizations. The red line marks the boundary of the topological phase.
	}
\end{figure}
Even with this high value of the disorder it is clear that the energy of the in-gap states is rather robust under this perturbation. For $t/\Delta=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 $\mu$ and in the energy of the first extended state. For higher ratios of $t/\Delta$, closer to the value of $N+1$ (cf. Eq.~\eqref{equation: criteria/ number of solutions for edge states even N} and the discussion under Eq.~\eqref{equation: spectrum: mu=0, even N, edges, cosh/cosh version}),
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 Fig.~\ref{figure: disorder-noaverage} 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 $\mu$. The existence of these crossings is in fact protected against local perturbations and they correspond to a switch of the fermionic parity \cite{Hegde-2015}.

\section{Conclusion}
\label{section: 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 ($\mu=0$, $t/\Delta = \pm 1$ in the notation of this paper), is easy to demonstrate. However, matters become soon complicated when generic values of the three parameters $\mu$, $\Delta$ 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 [\onlinecite{Kao-2014}], whereby for finite chemical potential ($\mu\neq0$) zero energy states only exists for discrete sets of $\mu (\Delta,\,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 $\xi$).

Beside of the Majorana lines, but still inside the topologial 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/\Delta$ , $\mu/\Delta$ and of the system size $L=Nd$. We found that with good accuracy the maximum energy $E_{0,max} \propto \exp(L/\xi)$. This energy, and hence the energy of all the boundary states, tends to zero in the thermodynamic limit $N\to \infty$. For fixed $N$  the ratio $t/\Delta$ can be varied until the decay length $\xi$ 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.

Although our treatment using Tetranacci polynomials for $\mu\neq 0$ is general (Secs. \ref{section: influence of non zero chemical potential}, \ref{section: MZM paths in the topological phase diagram for finite N}), we have dedicated special attention to two parameter choices in which the Tetranacci polynomials reduce to the generalized Fibonacci polynomials. The first case if that of zero chemical potential discussed in Sec. \ref{section: spectral analysis}, \ref{section: Eigenvectors and symmetries} 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 $\mu$ 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. At the Majorana lines the fact that $E=0$ decouples the two Majorana sublattices, and again allows us to use the simpler Fibonacci polynomials.

\section{Acknowledgments}

The authors thank the Elite Netzwerk Bayern for financial support via the IGK "Topological Insulators" and the Deutsche Forschungsgemeinschaft via SFB 1277 Project B04. We acknowledge useful discussions with A. Donarini, C. de Morais Smith and M. Wimmer.