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

The Kitaev chain is a one dimensional model based on a lattice of N spinless fermions. It is characterized by three parameters: the chemical potential μ, the hopping amplitude t, and the p-wave superconducting pairing constant Δ. The Kitaev Hamiltonian, written in a set of standard fermionic operators ${d}_{j},\;{d}_{j}^{{\dagger}}$, is [1, 2]

$$\begin{equation}{\hat{H}}_{\text{KC}}=-\mu \sum _{j=1}^{N}\;{d}_{j}^{{\dagger}}{d}_{j}+\sum _{j=1}^{N-1}\;\left({\Delta}\;{d}_{j}^{{\dagger}}{d}_{j+1}^{{\dagger}}\;-t\;{d}_{j+1}^{{\dagger}}{d}_{j}+\text{h.c.}\right),\end{equation} \tag{ 1 }$$

where the p-wave character allows interactions between particles of the same spin. The spin is thus not explicitly included in the following. We consider Δ and t to be real parameters from now on.

The Hamiltonian in equation (1) has drawn particular attention in the context of topological superconductivity, due to the possibility of hosting MZM at its end in a particular parameter range [2]. This can be seen by expressing the Kitaev Hamiltonian in terms of so called Majorana operators γ^A,B,

$$\begin{equation}\left(\begin{pmatrix}\hfill {d}_{j}\hfill \\ \hfill {d}_{j}^{{\dagger}}\hfill \end{pmatrix}\right)=:\;\frac{1}{\sqrt{2}}\;\left(\begin{pmatrix}\hfill 1\hfill & \hfill i\hfill \\ \hfill 1\hfill & \hfill -i\hfill \end{pmatrix}\right)\;\left(\begin{pmatrix}\hfill {\gamma }_{j}^{A}\hfill \\ \hfill {\gamma }_{j}^{B}\hfill \end{pmatrix}\right),\quad {\left({\gamma }^{A,B}\right)}^{{\dagger}}={\gamma }^{A,B},\end{equation} \tag{ 2 }$$

yielding the form

$$\begin{align}\hfill {\hat{H}}_{\text{KC}}& =-i\;\mu \;\sum _{j=1}^{N}\;{\gamma }_{j}^{A}{\gamma }_{j}^{B}+i\left({\Delta}+t\right)\sum _{j=1}^{N-1}{\gamma }_{j}^{B}{\gamma }_{j+1}^{A}\hfill \\ \hfill & \quad +\;i\left({\Delta}-t\right)\sum _{j=1}^{N-1}{\gamma }_{j}^{A}{\gamma }_{j+1}^{B}.\hfill \end{align} \tag{ 3 }$$

Notice that, in virtue of equation (2) it holds $\left\{{\gamma }_{j}^{A},\;{\left({\gamma }_{j}^{A}\right)}^{{\dagger}}\right\}=2\;{\left({\gamma }_{j}^{A}\right)}^{2}=1$, and similarly for ${\gamma }_{j}^{B}$. For the particular parameter settings Δ = ±t and μ = 0, which we call the Kitaev points, equation (3) leads to a 'missing' fermionic quasiparticle q_±:

$$\begin{equation}{q}_{+}=\frac{1}{\sqrt{2}}\left({\gamma }_{1}^{A}+i\;{\gamma }_{N}^{B}\right)\quad \left[{\Delta}=t\right],\end{equation} \tag{ 4a }$$

$$\begin{equation}{q}_{-}=\frac{1}{\sqrt{2}}\left({\gamma }_{1}^{B}+i\;{\gamma }_{N}^{A}\right)\quad \left[{\Delta}=-t\right].\end{equation} \tag{ 4b }$$

This quasiparticle has zero energy and is composed of two isolated Majorana states localised at the ends of the chain. In general, the condition of hosting MZM does not restrict to the Kitaev points (μ = 0, Δ = ±t). Further information on the existence of boundary modes is evinced from the bulk spectrum and the associated topological phase diagram. If the boundary modes have exactly zero energy, their Majorana nature can be proven by showing that they are eigenstates of the particle–hole operator $\mathcal{P}$. Equivalently, if ${\gamma }_{M}^{{\dagger}}$ is an operator associated to such an MZM it satisfies ${\gamma }_{M}^{{\dagger}}={\gamma }_{M}$ and ${\gamma }_{M}^{2}=1/2$.

The topological phase diagram is shortly reviewed in section 2.3.

The Hamiltonian from equation (1) in the limit of N → ∞ reads in k space

$$\begin{equation}{\hat{H}}_{\text{KC}}=\frac{1}{2}\sum _{k}\;{\hat{\psi }}_{k}^{{\dagger}}\;\mathcal{H}\left(k\right){\hat{\psi }}_{k},\quad {\hat{\psi }}_{k}={\left({d}_{k},\;{d}_{-k}^{{\dagger}}\right)}^{\text{T}},\end{equation} \tag{ 5 }$$

where we introduced the operators ${d}_{k}=\frac{1}{\sqrt{N}}\sum _{j}{\mathrm{e}}^{-\text{i}\;j\;kd}\;{d}_{j}$ and k lies inside the first Brillouin zone, i.e. $k\in \left[-\frac{\pi }{d},\frac{\pi }{d}\right]$ and d is the lattice constant. The 2 × 2 Bogoliubov–de Gennes (BdG) matrix

$$\begin{equation}\mathcal{H}\left(k\right)=\left[\begin{bmatrix}{cc}\hfill -\mu -2t\;\mathrm{cos}\left(kd\right)\hfill & \hfill -2i{\Delta}\mathrm{sin}\left(kd\right)\hfill \\ \hfill 2i{\Delta}\;\mathrm{sin}\left(kd\right)\hfill & \hfill \mu +2t\;\mathrm{cos}\left(kd\right)\hfill \end{bmatrix}\right]\end{equation} \tag{ 6 }$$

is easily diagonalized thus yielding the excitation spectrum

$$\begin{equation}{E}_{{\pm}}\left(k\right)={\pm}\sqrt{4{{\Delta}}^{2}\;{\mathrm{sin}}^{2}\left(kd\right)+{\left[\mu +2t\;\mathrm{cos}\left(kd\right)\right]}^{2}}.\end{equation} \tag{ 7 }$$

Note that for μ = 0 equation (7) predicts a gapped spectrum whose gap width is either 4Δ (|Δ| < |t|) or 4t (|t| < |Δ|).

The BdG Hamiltonian (6) 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}=\mathbb{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 [28]. Note that all symmetries square to +1, placing the Kitaev chain in the BDI class [31].

The winding number is given by [28, 32]

$$\begin{equation}\nu =\frac{1}{2\pi }\underset{-\pi /d}{\overset{\pi /d}{\int }}\mathrm{d}k\;{\partial }_{k}\;w\left(k\right),\end{equation} \tag{ 8 }$$

where $w\left(k\right)=\mathrm{arg}\left[2{\Delta}\;\mathrm{sin}\left(kd\right)+i\;\left(\mu +2t\;\mathrm{cos}\left(kd\right)\right)\right]$ and ∂_kw(k) is the winding number density. A trivial phase corresponds to ν = 0, a non trivial one to finite integer values of ν. The winding number relates bulk properties to the existence of boundary (not necessarily MZM) states in a finite chain. This property is known as bulk-edge correspondence. Due to their topological nature, their existence is robust against small perturbations, like disorder. This point is further discussed in section 7.

The phase diagram constructed using the winding number invariant is shown in figure 1. The meaning of two different values for the winding number is clearer when we recall the Kitaev Hamiltonian in the Majorana basis. In a finite chain the leftmost lattice site consists of the A Majorana operator ${\gamma }_{1}^{A}$ connected to the bulk by the i(Δ − t) hopping and the B Majorana operator ${\gamma }_{1}^{B}$ connected by the i(Δ + t) hopping. With Δ > 0 and t < 0 (the ν = +1 phase) the Majorana state at the left end of the chain will consist mostly of the weakly connected ${\gamma }_{1}^{B}$. If t > 0 (the ν = −1 phase), ${\gamma }_{1}^{A}$ is connected to the bulk more weakly and contributes most to the left end bound state.

Zoom In Zoom Out Reset image size

The boundaries between different topological phases can be obtained from the condition of closing the bulk gap, i.e. E_±(k) = 0 (cf equation (7)). That is only possible if both terms under the square root vanish. The condition of Δ ≠ 0 forces the gap closing to occur at kd = 0 or k = πd, and the remaining term vanishes at these momenta if μ = ±2t. The four insets in figure 1 show the behavior of w(k), leading to either a zero (for w(−π) = w(π)) or non zero winding number, see equation (8).

Physically speaking, the Kitaev chain is in the topological phase provided that Δ ≠ 0 and the chemical potential lies inside the 'normal' band (|μ| ⩽ 2|t|).

The BdG Hamiltonian is block diagonal and its characteristic polynomial ${P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)=\mathrm{det}\left[\lambda \;\mathbb{1}-{\mathcal{H}}_{\text{KC}}\right]$ factorises as

$$\begin{equation}{P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)={P}_{\lambda }\left(C\right)\;{P}_{\lambda }\left(-C\right).\end{equation} \tag{ 13 }$$

The tridiagonal structure of C straightforwardly yields the spectrum of a normal conducting, linear chain [33]

$$\begin{equation}{E}_{{\pm}}^{{\Delta}=0}\left({k}_{n}\right)={\pm}\left[\mu +2t\;\mathrm{cos}\left({k}_{n}d\right)\right],\quad {k}_{n}d=\frac{n\pi }{N+1},\end{equation} \tag{ 14 }$$

where n runs from 1 to N. Since ${k}_{n}\in \mathbb{R}$, only bulk states exist for Δ = 0.

In the beginning we consider both t and μ to be zero and include μ ≠ 0 in a second step. The parameter setting leads to a vanishing matrix C, see equation (11), and the characteristic polynomial of the system reads:

$$\begin{equation}{P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)=\mathrm{det}\left[\begin{bmatrix}{cc}\hfill \lambda \;\mathbb{1}\hfill & \hfill -S\hfill \\ \hfill S\hfill & \hfill \lambda \;\mathbb{1}\hfill \end{bmatrix}\right],\end{equation} \tag{ 15 }$$

where we used the property S^† = −S. Due to the fact that the commutator $\left[\mathbb{1},\;S\right]=0$ vanishes⁴, one finds [34]

$$\begin{equation}{P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)=\mathrm{det}\left({\lambda }^{2}\;\mathbb{1}+{S}^{2}\right).\end{equation} \tag{ 16 }$$

The characteristic polynomial can still be simplified to the product

$$\begin{equation}{P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)={P}_{\lambda }\left(iS\right)\;{P}_{\lambda }\left(-iS\right).\end{equation} \tag{ 17 }$$

The matrix iS is hermitian and describes a linear chain with hopping iΔ. As a consequence, we find the spectrum to be [33]

$$\begin{equation}{E}_{{\pm}}\left({k}_{n}\right)={\pm}\left[2{\Delta}\;\mathrm{cos}\left({k}_{n}d\right)\right],\quad {k}_{n}d=\frac{n\pi }{N+1},\end{equation} \tag{ 18 }$$

where n runs from 1 to N and each eigenvalue is twice degenerated. Notice the phase shift by π/2 compared to the spectrum of an infinite chain equation (7). We discuss this phase shift in more detail in section 3.3.

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

The chemical potential μ can be included easily. Exploiting the properties of ${\mathcal{H}}_{\text{KC}}$, we find the characteristic polynomial to be

$$\begin{align}\hfill {P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)& =\mathrm{det}\left[\begin{bmatrix}{cc}\hfill \left(\lambda +\mu \right)\;\mathbb{1}\hfill & \hfill -S\hfill \\ \hfill S\hfill & \hfill \left(\lambda -\mu \right)\;\mathbb{1}\hfill \end{bmatrix}\right]\hfill \\ \hfill & =\;\mathrm{det}\left[\left({\lambda }^{2}-{\mu }^{2}\right)\;\mathbb{1}+{S}^{2}\right]\hfill \\ \hfill & =\;\mathrm{det}\left[{{\Lambda}}^{2}\;\mathbb{1}+{S}^{2}\right],\hfill \end{align} \tag{ 19 }$$

with Λ² := λ² − μ². The same treatment as in the previous μ = 0 case yields ${P}_{\lambda }\left({\mathcal{H}}_{\text{KC}}\right)={P}_{{\Lambda}}\left(iS\right)\;{P}_{{\Lambda}}\left(-iS\right)$. Consequently the spectrum is

$$\begin{equation}{E}_{{\pm}}^{t=0}\left({k}_{n}\right)={\pm}\sqrt{{\mu }^{2}+4{{\Delta}}^{2}\;{\mathrm{cos}}^{2}\left({k}_{n}d\right)},\quad {k}_{n}d=\frac{n\pi }{N+1},\end{equation} \tag{ 20 }$$

where n runs again from 1 to N. Again no boundary modes are found for t = 0.

The calculation of the spectrum for μ = 0 requires a more technical approach, since the structure of the BdG Hamiltonian equation (10) prohibits standard methods.

One important feature of the Kitaev chain can be appreciated inspecting equation (3). The entire model is equivalent to two coupled SSH-like chains [30, 35] containing both the hopping parameters $a\;{:=}\;i\left({\Delta}-t\right)$ and $b\;{:=}\;i\left({\Delta}+t\right)$, see figure 2. Explicitly,

$$\begin{align}\hfill {\hat{H}}_{\text{KC}}& =\frac{1}{2}\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)+\text{h.c.}\hfill \\ \hfill & \quad +\;\frac{1}{2}\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)+\text{h.c.}\hfill \\ \hfill & \quad -\;i\mu \sum _{j=1}^{N}\;{\gamma }_{j}^{A}{\gamma }_{j}^{B},\hfill \end{align} \tag{ 21 }$$

where N_1,2 depend on N. If N is even we have N₁ = N/2 and N₂ = N₁ − 1, while N₁ = N₂ = (N − 1)/2 for odd N. Independent of the number of atoms, the first and the second lines in equation (21) describe two SSH-like chains, coupled by the chemical potential μ. We define here the SSH-like basis of the Kitaev chain as:

$$\begin{align}\hfill {\hat{{\Psi}}}_{\text{SSH}}^{\text{even}}& ={\left({\gamma }_{1}^{A},\;{\gamma }_{2}^{B},\dots ,\;{\gamma }_{N-1}^{A},\;{\gamma }_{N}^{B}\vert {\gamma }_{1}^{B},\;{\gamma }_{2}^{A},\dots ,\;{\gamma }_{N-1}^{B},\;{\gamma }_{N}^{A}\right)}^{T},\hfill \\ \hfill \\ \hfill {\hat{{\Psi}}}_{\text{SSH}}^{\text{odd}}& ={\left({\gamma }_{1}^{\text{A}},\;{\gamma }_{2}^{B},\dots ,\;{\gamma }_{N-1}^{B},\;{\gamma }_{N}^{A}\vert {\gamma }_{1}^{B},\;{\gamma }_{2}^{A},\dots ,\;{\gamma }_{N-1}^{A},\;{\gamma }_{N}^{B}\right)}^{\text{T}},\hfill \end{align} \tag{ 22 }$$

where '|' marks the boundary between both chains. We call the first one, starting always with ${\gamma }_{1}^{A}$, the α chain, and the second the β chain, such that ${{\Psi}}_{\text{SSH}}^{\text{even,odd}}={\left({\to {\gamma }}_{\alpha }\;\vert \;{\to {\gamma }}_{\beta }\right)}^{\text{T}}$. The BdG Hamiltonian in the SSH-like basis reads

$$\begin{equation}{\mathcal{H}}_{\text{KC}}^{\text{SSH}}=\left[\begin{bmatrix}{cc}\hfill {\mathcal{H}}_{\alpha }\hfill & \hfill \tau \hfill \\ \hfill {\tau }^{{\dagger}}\hfill & \hfill {\mathcal{H}}_{\beta }\hfill \end{bmatrix}\right],\end{equation} \tag{ 23 }$$

with ${\hat{H}}_{\text{KC}}=\frac{1}{2}{\hat{{\Psi}}}_{\text{SSH}}^{{\dagger}}{\mathcal{H}}_{\text{KC}}^{\text{SSH}}{\hat{{\Psi}}}_{\text{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 τ and τ^†, which contain only the chemical potential μ, in a diagonal arrangement specified below.

Zoom In Zoom Out Reset image size

The pattern of these matrices is slightly different for even and odd number of sites. If N is even we find

$$\begin{equation}{\mathcal{H}}_{\alpha }^{\text{even}}=\left[\begin{bmatrix}{cc}\hfill 0\hfill & \hfill a\hfill \\ \hfill -a\hfill & \hfill 0\hfill & \hfill b\hfill \\ \hfill \hfill & \hfill -b\hfill & \hfill 0\hfill & \hfill a\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -a\hfill & \hfill 0\hfill & \hfill b\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -b\hfill & \hfill 0\hfill & \hfill a\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -a\hfill & \hfill 0\hfill \end{bmatrix}\right],\end{equation} \tag{ 24 }$$

$$\begin{equation}{\mathcal{H}}_{\beta }^{\text{even}}=\left[\begin{bmatrix}{cc}\hfill 0\hfill & \hfill b\hfill \\ \hfill -b\hfill & \hfill 0\hfill & \hfill a\hfill \\ \hfill \hfill & \hfill -a\hfill & \hfill 0\hfill & \hfill b\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -b\hfill & \hfill 0\hfill & \hfill a\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -a\hfill & \hfill 0\hfill & \hfill b\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -b\hfill & \hfill 0\hfill \end{bmatrix}\right],\end{equation} \tag{ 25 }$$

and ${\tau }^{\text{even}}=-i\mu \;{\mathbb{1}}_{N/2}\otimes {\tau }_{z}$, where τ_z denotes the Pauli matrix. The odd N expressions are achieved by removing the last line and column in ${\mathcal{H}}_{\alpha }^{\text{even}}$, ${\mathcal{H}}_{\beta }^{\text{even}}$ and τ^even.

As shown in more detail in appendix B, for μ = 0 the characteristic polynomial can be expressed as the product of two polynomials of order N

$$\begin{equation}{P}_{\lambda }{\left({\mathcal{H}}_{\text{KC}}\right)}_{\mu =0}={\zeta }_{N}\left(\lambda ,\;a,\;b\right)\;\;{{\epsilon}}_{N}\left(\lambda ,\;a,\;b\right),\end{equation} \tag{ 26 }$$

where the product form reflects the fact that the Kitaev chain is given in terms of two uncoupled SSH-like chains, as illustrated in figure 2. Even though the polynomials ζ_N and _N belong to different SSH-like chains, both obey a common recursion formula typical of Fibonacci polynomials [36–38]

$$\begin{equation}{\zeta }_{j+2}=\left[{\lambda }^{2}+{a}^{2}+{b}^{2}\right]\;{\zeta }_{j}-{a}^{2}{b}^{2}\;{\zeta }_{j-2},\end{equation} \tag{ 27 }$$

and differ only in their initial values

$$\begin{equation}\left(\begin{pmatrix}\hfill {\zeta }_{-1}\hfill \\ \hfill {\zeta }_{0}\hfill \\ \hfill {\zeta }_{1}\hfill \\ \hfill {\zeta }_{2}\hfill \end{pmatrix}\right)=\left(\begin{pmatrix}0\hfill \\ 1\hfill \\ \lambda \hfill \\ \hfill {\lambda }^{2}+{b}^{2}\hfill \end{pmatrix}\right),\;\left(\begin{pmatrix}\hfill {{\epsilon}}_{-1}\hfill \\ {{\epsilon}}_{0}\hfill \\ {{\epsilon}}_{1}\hfill \\ {{\epsilon}}_{2}\hfill \end{pmatrix}\right)=\left(\begin{pmatrix}0\hfill \\ 1\hfill \\ \lambda \hfill \\ \hfill {\lambda }^{2}+{a}^{2}\hfill \end{pmatrix}\right).\end{equation} \tag{ 28 }$$

Fundamental properties of Fibonacci polynomials are summarized in appendix A. The common sublattice structure of both chains sets the stage for a relationship between ζ_j and _j: the exchange of a's and b's enables us to pass from one to the other

$$\begin{equation}{\zeta }_{j}\left(\lambda ,\;a,\;b\right)={{\epsilon}}_{j}\left(\lambda ,\;b,\;a\right),\quad \forall \;j.\end{equation} \tag{ 29 }$$

Moreover, equation (27) implies that a Kitaev chain with even number of sites N is fundamentally different from the one with an odd number of sites. This property is a known feature of SSH chains [39]. The difference emerges since, according to equations (B20), (B22), it holds

$$\begin{equation}{\zeta }_{\text{odd}}\left(\lambda ,\;a,\;b\right)={{\epsilon}}_{\text{odd}}\left(\lambda ,\;a,\;b\right),\end{equation} \tag{ 30 }$$

because the number of a and b type bondings in both subchains is the same. This leads to twice degenerate eigenvalues. An equivalent relationship for even N does not exist. The closed form for ζ_j and _j, as well as their factorization, is derived in appendix B.

The characteristic polynomial can be used to obtain the determinant of the Kitaev chain, here for μ = 0, because evaluating it at λ = 0 leads to:

$$\begin{equation*}{P}_{\lambda =0}{\left({\mathcal{H}}_{\text{KC}}\right)}_{\mu =0}=\mathrm{det}{\left({\mathcal{H}}_{\text{KC}}\right)}_{\mu =0}.\end{equation*}$$

According to equation (26) we need only to know ζ_N and _N at λ = 0. The closed form expression for ζ_j at λ = 0 reduces to

$$\begin{equation}{\zeta }_{j}{\vert }_{\lambda =0}=\left\{\begin{matrix}\hfill 0,\hfill & \hfill \text{if}\;\;j\;\;\;\text{is}\;\text{odd}\hfill \\ \hfill {b}^{\;j},\hfill & \hfill \;\text{else}\;\hfill \end{matrix}\right.,\end{equation} \tag{ 31 }$$

while _j|_λ=0. follows from equation (29). We find that there are always zero energy eigenvalues for odd N, but not in general for even N, as it follows from

$$\begin{equation}\mathrm{det}\left({\mathcal{H}}_{\text{KC}}^{\mu =0}\right)=\left\{\begin{matrix}\hfill 0,\hfill & \hfill N\;\text{odd}\hfill \\ \hfill {\left[{{\Delta}}^{2}-{t}^{2}\right]}^{N},\hfill & \hfill N\;\text{even}\hfill \end{matrix}\right..\end{equation} \tag{ 32 }$$

Additional features of the spectrum are discussed in the following.

3.3.1. Odd N

The spectrum for odd N is given by two contributions

$$\begin{equation}{E}_{{\pm}}^{\mu =0}=0,\quad \left(\text{twofold}\right)\end{equation} \tag{ 33 }$$

$$\begin{equation}{E}_{{\pm}}^{\mu =0}\left({k}_{n}\right)={\pm}\sqrt{4{{\Delta}}^{2}\;{\mathrm{sin}}^{2}\left({k}_{n}d\right)+4{t}^{2}\;{\mathrm{cos}}^{2}\left({k}_{n}d\right)},\end{equation} \tag{ 34 }$$

where k_nd = nπ/(N + 1) and n runs from 1 to N, except for n = (N + 1)/2. This constraint on n is a consequence of the equations (67), (69) below which show that the boundary condition equation (70) cannot be satisfied for kd = π/2. Hence, no standing wave can be formed.

Each zero eigenvalue belongs to one chain. As discussed below, two decaying states are associated to equation (33), whose wave functions are discussed in section 4.2. These states are MZM.

3.3.2. Even N

In the situation of even N we find for the Kitaev's bulk spectrum at zero μ

$$\begin{equation}{E}_{{\pm}}^{\mu =0}\left(k\right)={\pm}\sqrt{4{{\Delta}}^{2}\;{\mathrm{sin}}^{2}\left(kd\right)+4{t}^{2}\;{\mathrm{cos}}^{2}\left(kd\right)},\end{equation} \tag{ 35 }$$

where the momenta k are in general not equidistant in the first Brillouin zone. Rather, the bulk quantization condition follows from the interplay between Δ and t and is captured in form of the functions f_β,α(k) (cf appendix B.2),

$$\begin{equation}{f}_{\beta ,\alpha }\left(k\right){:=}\mathrm{tan}\left[kd\left(N+1\right)\right]\;{\pm}\frac{{\Delta}}{t}\;\mathrm{tan}\left(kd\right),\end{equation} \tag{ 36 }$$

whose zeros

$$\begin{equation}{f}_{\beta ,\alpha }\left(k\right)\;!{=}\;0,\quad kd\ne 0,\;\pi /2\end{equation} \tag{ 37 }$$

define the allowed values of k. Note that kd = 0, π/2 are excluded as solutions, due to their trivial character. The functions f_β,α(k) follow from the factorisation of the polynomials _N and ζ_N. The negative sign in equation (36) belongs to the α subchain, while the positive one to the β subchain. The spectrum following from equation (37) is illustrated in figure 3. We observe that equations (35) and (37) hold for all values of t and Δ, independent of whether |Δ| is larger or smaller than |t|. The two situations are connected by a phase shift of the momentum kd → kd + π/2, which influences both the spectrum and the quantization condition. In the end all different ratios of Δ and t are captured by equations (35) and (37), due to the periodicity of the spectrum.

Zoom In Zoom Out Reset image size

However, when we consider decaying or edge states this periodicity is lost (see equations (40) and (41) below) and |t| ≶ |Δ| lead to different quantization rules. The hermiticity of the Hamiltonian allows a pure imaginary momentum for μ = 0, but a simple exchange of k to iq in equation (36) does not lead to the correct results. We introduce here the functions

$$\begin{equation}{h}_{\beta ,\alpha }\left(q\right){:=}\mathrm{tanh}\left[qd\left(N+1\right)\right]\;{\pm}\;m\;\mathrm{tanh}\left(qd\right),\end{equation} \tag{ 38 }$$

similar to f_β,α(k) in equation (36), where m contains both ratios of Δ and t:

$$\begin{equation}m\;{:=}\;\left\{\begin{matrix}\hfill \frac{{\Delta}}{t},\quad \;\mathrm{i}\mathrm{f}\;\vert {\Delta}\vert {\geqslant}\vert t\vert \hfill \\ \hfill \frac{t}{{\Delta}},\quad \;\mathrm{i}\mathrm{f}\;\vert t\vert {\geqslant}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation} \tag{ 39 }$$

Again, the positive sign in equation (38) belongs to the β chain and the negative one to the α chain. The exact quantization criterion is provided by the zeros of h_β,α(q),

$$\begin{equation}{h}_{\beta ,\alpha }\left(q\right)\;!{=}0,\quad q\ne 0,\end{equation} \tag{ 40 }$$

as illustrated in figure 4. The associated energies follow from the dispersion relation

$$\begin{equation}E\left(q\right)={\pm}\;\left\{\begin{matrix}\hfill \sqrt{4{t}^{2}\;{\mathrm{cosh}}^{2}\left(qd\right)-4{{\Delta}}^{2}\;{\mathrm{sinh}}^{2}\left(qd\right)},\quad \mathrm{i}\mathrm{f}\;\vert {\Delta}\vert {\geqslant}\vert t\vert \hfill \\ \hfill \sqrt{4{{\Delta}}^{2}\;{\mathrm{cosh}}^{2}\left(qd\right)-4{t}^{2}\;{\mathrm{sinh}}^{2}\left(qd\right)},\quad \mathrm{i}\mathrm{f}\;\vert t\vert {\geqslant}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation} \tag{ 41 }$$

We notice that equation (41) is only well defined for zero or positive arguments of the square root. Indeed, all solutions of equation (40), if existent, lie always inside this range, because using h_β,α(q) = 0 in equation (41) yields

$$\begin{equation}E\left(q\right)={\pm}\;2\;\frac{\mathrm{cosh}\left(qd\right)}{\mathrm{cosh}\left[qd\left(N+1\right)\right]}\;\cdot \;\mathrm{min}\left\{\vert {\Delta}\vert ,\vert t\vert \right\}.\end{equation} \tag{ 42 }$$

Hence, each wavevector from equation (40) 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 equation (42) is always smaller than one. We can restrict ourselves to find only positive solutions qd, due to the time reversal symmetry. The number of physically different solutions of equation (40) is zero or two and it follows always from the equation containing the positive factor m or −m. Consequently, according to equation (38), only none or two gap modes can form and both belong to the same subchain, α or β. Moreover a solution exists if, and only if, $\vert m\vert \in \left[1,\;N+1\right]$.

Zoom In Zoom Out Reset image size

In the limiting case when |m| → 1, i.e. at the Kitaev points, the solution qd → ∞ and the associated energies E_± from equation (42) go to zero. The eigenstate will be a Majorana zero energy mode, see section 4.1.

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

Increasing the value of |m| beyond N + 1 entails the absence of imaginary solutions. The number of eigenvalues of a Kitaev chain is still 2N for a fixed number of sites and consequently equation (37) leads now to N real values for kd, instead of N − 1. In other words, the two former gap modes have moved to two extended states and their energy lies now within the bulk region of the spectrum, even though the system is still fully gaped. This effect holds for the Kitaev chain as well as for SSH chains. Physically this means, that a 'boundary' mode with imaginary momentum q and corresponding decay length ξ ∝ 1/q reached the highest possible delocalisation in the chain.

The limit of N → ∞ yields always two zero energy boundary modes; since the momentum is qd = arctanh(|1/m|), due to equations (38), (40) and according to equation (42) the energy goes to zero. If we consider the odd N situation in the limit of an infinite number of sites, we have there two zero energy boundary modes as well. The results of this section are summarized in table 1.

Table 1. Overview of the quantization rule for the wave vectors of the finite Kitaev chain in different scenarios. The wavevectors k, k_n, κ_1,2 used together with equation (7) and q with equation (41) yield the correct finite system energies and $k,\;{k}_{n},q\in \mathbb{R}$, ${\kappa }_{1,2}\in \mathbb{C}$. Notice: n = 1, ..., N and $m=\frac{t}{{\Delta}}$ ($m=\frac{{\Delta}}{t}$) for |t| ⩾ |Δ| (|Δ| > |t|).

Requirements	Quantisation rule	Zero modes	Equation for eigenstate elements	Majorana character
Δ = 0:	${k}_{n}d=\frac{n\pi }{N+1}$	Yes, if for some n:		No
			$\mu ={\mu }_{n}=2t\;\mathrm{cos}\left(\frac{n\pi }{N+1}\right)$
t = 0:	${k}_{n}d=\frac{\pi }{2}+\frac{n\pi }{N+1}$	No		No
μ = 0, N odd:	${k}_{n}d=\frac{n\pi }{N+1},n\ne \frac{N+1}{2}$	No	(67)–(69)	No
	$qd=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(1/\vert m\vert \right)$	Yes	(71), (72)	Yes
			(75)–(78)	Yes
μ = 0, N even:	$\mathrm{tan}\left[kd\left(N+1\right)\right]=\mp \frac{{\Delta}}{t}\mathrm{tan}\left(kd\right)$	No	(51), (52)	No
	$\mathrm{tanh}\left[qd\left(N+1\right)\right]=\mp m\;\mathrm{tanh}\left(qd\right)$	Only if Δ = ±t	(61)–(64)	Yes
		Otherwise	(51), (52)	No
$t,\;{\Delta},\;\mu \in \mathbb{R}$	$\frac{{\mathrm{sin}}^{2}\left[\frac{{\kappa }_{1}+{\kappa }_{2}}{2}\left(N+1\right)\right]}{{\mathrm{sin}}^{2}\left[\frac{{\kappa }_{1}-{\kappa }_{2}}{2}\left(N+1\right)\right]}=\frac{1+{\left(\frac{{\Delta}}{t}\right)}^{2}\;{\mathrm{cot}}^{2}\left(\frac{{\kappa }_{1}-{\kappa }_{2}}{2}\right)}{1+{\left(\frac{{\Delta}}{t}\right)}^{2}\;{\mathrm{cot}}^{2}\left(\frac{{\kappa }_{1}+{\kappa }_{2}}{2}\right)}$	only for $\mu ={\mu }_{n}\in \mathbb{R}$ and	(100), (102)	No
		${\mu }_{n}=2\sqrt{{t}^{2}-{{\Delta}}^{2}}\mathrm{cos}\left(\frac{n\pi }{N+1}\right)$	(109), (110)	Yes
			(112), (113)	Yes

We define the vectors ${\to {v}}_{\alpha }$ and ${\to {v}}_{\beta }$ via the entries

$$\begin{equation}{\to {v}}_{\alpha }={\left({x}_{1},\;{y}_{1},\;{x}_{2},\;{y}_{2},\dots ,\;{x}_{N/2},\;{y}_{N/2}\right)}^{\text{T}},\end{equation} \tag{ 45 }$$

$$\begin{equation}{\to {v}}_{\beta }={\left({\mathcal{X}}_{1},\;{\mathcal{Y}}_{1},\;{\mathcal{X}}_{2},\;{\mathcal{Y}}_{2},\dots ,\;{\mathcal{X}}_{N/2},\;{\mathcal{Y}}_{N/2}\right)}^{\text{T}},\end{equation} \tag{ 46 }$$

where x, y and $\mathcal{X},\;\mathcal{Y}$ are associated to the A and B sublattices, respectively. The internal structure of ${\to {v}}_{\alpha }$ (${\to {v}}_{\beta }$) reflects the unit cell of an SSH-like chain and thus simplifies the calculation. In the real space x_l (${\mathcal{X}}_{l}$) belongs to site j = 2l − 1 and y_l (${\mathcal{Y}}_{l}$) to j = 2l, where j = 1, ..., N.

Searching for solutions on the subchain α implies setting ${\to {v}}_{\beta }=\to {0}$ and solving $\left({\mathcal{H}}_{\alpha }^{\text{even}}-{E}_{{\pm}}{\mathbb{1}}_{N}\right){\to {v}}_{\alpha }=\to {0}$. The elements of ${\to {v}}_{\alpha }$ obey

$$\begin{equation}\;a\;{y}_{1}={E}_{{\pm}}\;{x}_{1},\end{equation} \tag{ 47 }$$

$$\begin{equation}-a\;{x}_{N/2}={E}_{{\pm}}\;{y}_{N/2},\end{equation} \tag{ 48 }$$

and

$$\begin{equation}b\;{x}_{l+1}-a\;{x}_{l}={E}_{{\pm}}\;{y}_{l},\end{equation} \tag{ 49 }$$

$$\begin{equation}a\;{y}_{l+1}-b\;{y}_{l}={E}_{{\pm}}\;{x}_{l+1},\end{equation} \tag{ 50 }$$

where l runs from 1 to N/2 − 1. The solution for Δ ≠ ±t is (in agreement with [40])

$$\begin{equation}\frac{{y}_{l}}{{x}_{1}}=\frac{{E}_{{\pm}}\left(\theta \right)}{a}\;{T}_{l}^{\text{e}}\left(\theta \right),\end{equation} \tag{ 51 }$$

$$\begin{equation}\frac{{x}_{l}}{{x}_{1}}={T}_{l}^{\text{e}}\left(\theta \right)-\frac{b}{a}{T}_{l-1}^{\text{e}}\left(\theta \right),\end{equation} \tag{ 52 }$$

where l = 1, ..., N/2, θ/(2d) denotes the momentum k (q) for extended (gap) states and E_± is the dispersion relation associated to k (equation (35)), or q (equation (41)). The entries of the eigenvectors are essentially sine functions for the extended states

$$\begin{equation}{T}_{l}^{\text{e}}\left(k\right){:=}\;\frac{\mathrm{sin}\left(2\;kd\;l\right)}{\mathrm{sin}\left(2\;kd\right)},\end{equation} \tag{ 53 }$$

and hyperbolic sine functions for the decaying states

$$\begin{equation}{T}_{l}^{\text{e}}\left(q\right){:=}{s}^{l-1}\;\;\;\frac{\mathrm{sinh}\left(2\;qd\;l\right)}{\mathrm{sinh}\left(2\;qd\right)},\end{equation} \tag{ 54 }$$

where the prefactor s depends on the ratio of Δ and t:

$$\begin{equation*}s=\left\{\begin{matrix}\hfill +1,\hfill & \hfill \vert {\Delta}\vert { >}\vert t\vert \hfill \\ \hfill -1,\hfill & \hfill \vert t\vert { >}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation*}$$

An illustration of ${\to {\psi }}^{\alpha }$ is given in figure 5. The allowed momenta k or q follow from the open boundary conditions

$$\begin{equation}{y}_{0}={x}_{\frac{N}{2}+1}=0.\end{equation} \tag{ 55 }$$

The first condition is satisfied due to T₀(θ) = 0 for any momentum. The second condition yields the quantization rules f_α(k) = 0 and h_α(q) = 0 for the α chain, see equations (37) and (40).

The eigenvector ${\to {\psi }}^{\beta }$ entails ${\to {v}}_{\alpha }=0$ and the entries of ${\to {v}}_{\beta }$ follow essentially by replacing a's and b's in the equations (51) and (52). We find

$$\begin{equation}\frac{{\mathcal{Y}}_{l}}{{\mathcal{X}}_{1}}=\frac{{E}_{{\pm}}}{b}\;{T}_{l}^{\text{e}}\left(\theta \right),\end{equation} \tag{ 56 }$$

$$\begin{equation}\frac{{\mathcal{X}}_{l}}{{\mathcal{X}}_{1}}={T}_{l}^{\text{e}}\left(\theta \right)-\frac{a}{b}{T}_{l-1}^{\text{e}}\left(\theta \right),\end{equation} \tag{ 57 }$$

where j = 1, ..., N/2 and Δ ≠ ±t. The quantisation condition follows from the open boundary condition:

$$\begin{equation*}{\mathcal{Y}}_{0}=0,\quad {\mathcal{X}}_{\frac{N}{2}+1}=0,\end{equation*}$$

and k (q) obey f_β(k) = 0 (h_β(q) = 0). Further, from the quantization rules it follows that gap modes belong always to the same subchain α or β for even N.

Zoom In Zoom Out Reset image size

As illustrated in figures 5–7 our states are symmetric w.r.t. the center of the SSH-like chains. This symmetry is visible in alternative versions of the equations (52) and (57) (l = 1, ..., N/2), whereby ${\mathit{x}}_{\frac{N}{2}+1-l}={T}_{l}^{\text{e}}\left(\theta \right)\;{\mathit{x}}_{\frac{N}{2}}$, which holds for all eigenstates. Together with equations (51), we find in general

$$\begin{equation}{x}_{\frac{N}{2}+1-l}={y}_{l}\;\frac{{x}_{N/2}}{{y}_{1}}\end{equation} \tag{ 58 }$$

and similarly ${\mathcal{X}}_{\frac{N}{2}+1-l}={\mathcal{Y}}_{l}\;{\mathcal{X}}_{\frac{N}{2}}/{\mathcal{Y}}_{1}$. Recalling the definition of the SSH-like basis, equation (22), and introducing the operators ${\psi }_{\alpha ,\beta }^{{\dagger}}$ associated to the states ${\to {\psi }}^{\alpha ,\beta }$ in equation (44), we find the expression

$$\begin{equation}{\psi }_{\alpha }^{{\dagger}}=\frac{1}{{v}_{\alpha }}\left[\sum _{j=1}^{N/2}\;{x}_{j}\;{\left({\gamma }_{2j-1}^{A}\right)}^{{\dagger}}+\sum _{j=1}^{N/2}\;{y}_{j}\;{\left({\gamma }_{2j}^{B}\right)}^{{\dagger}}\right],\end{equation} \tag{ 59 }$$

where v_α is the norm of the vector ${\to {v}}_{\alpha }$. A similar term is found for ${\psi }_{\beta }^{{\dagger}}$. We notice that equations (59) and (60) below are true for all kinds of eigenstates, i.e. extended, decaying states and MZM, of the BdG Hamiltonian in equation (23) at μ = 0. The character (statistics) of the operators depends on whether ${\left({\psi }_{\alpha ,\beta }^{{\dagger}}\right)}^{2}$ is 0 or 1/2. The property $\left\{{\gamma }_{j}^{r},{\gamma }_{k}^{s}\right\}={\delta }_{j,k}{\delta }_{r,s}$ with (r, s ∈ {A, B}) yields

$$\begin{equation}{\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}=\frac{1}{2\;{v}_{\alpha }^{2}}\sum _{j=1}^{N/2}\;\left({x}_{j}^{2}+{y}_{j}^{2}\right).\end{equation} \tag{ 60 }$$

The symmetry in equations (58) states that ${\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}$ is essentially determined by ${\left({x}_{N/2}/{y}_{1}\right)}^{2}$. For a ≠ 0 and consequently E ≠ 0, we find from equations (47), (48) and (58) that ${\left({x}_{N/2}/{y}_{1}\right)}^{2}=-1$, which yields ${\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}=0$. Thus the operators associated to the finite energy states ${\to {\psi }}^{\alpha }$, including the ones depicted in figures 5(a) and 6(a), obey fermionic statistics. This result holds also true in the case a = 0 and E ≠ 0 as can be seen by using the corresponding eigenstates (appendix C). Similar results hold for ${\left({\psi }_{\beta }^{{\dagger}}\right)}^{2}$.

Zoom In Zoom Out Reset image size

We turn now to Majorana zero modes, which at μ = 0 only exist at the Kitaev points Δ = ±t, see equation (32).

When Δ = t we find two zero energy modes ${\to {\psi }}_{A}^{\alpha }=\left(\begin{pmatrix}\hfill {\to {v}}_{\alpha ,A}\hfill \\ \to {0}\hfill \end{pmatrix}\right)$, ${\to {\psi }}_{B}^{\alpha }=\left(\begin{pmatrix}\hfill {\to {v}}_{\alpha ,B}\hfill \\ \to {0}\hfill \end{pmatrix}\right)$ each localised at one end of the α chain:

$$\begin{equation}{\to {v}}_{\alpha ,A}={\left(1,\;0,\;0,\dots ,\;0\right)}^{\text{T}},\end{equation} \tag{ 61 }$$

$$\begin{equation}{\to {v}}_{\alpha ,B}={\left(0,\;0,\dots ,\;0,\;1\right)}^{\text{T}},\end{equation} \tag{ 62 }$$

and ${\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}=1/2$ in equation (60). In contrast, both zero energy modes are on the β chain for Δ = −t. We find ${\to {\psi }}_{A}^{\beta }=\left(\begin{pmatrix}\to {0}\hfill \\ \hfill {\to {v}}_{\beta ,A}\hfill \end{pmatrix}\right)$, ${\to {\psi }}_{B}^{\beta }=\left(\begin{pmatrix}\to {0}\hfill \\ \hfill {\to {v}}_{\beta ,B}\hfill \end{pmatrix}\right)$ with

$$\begin{equation}{\to {v}}_{\beta ,\;B}={\left(1,\;0,\;0,\dots ,\;0\right)}^{\text{T}},\end{equation} \tag{ 63 }$$

$$\begin{equation}{\to {v}}_{\beta ,\;A}={\left(0,\;0,\dots ,\;0,\;1\right)}^{\text{T}}.\end{equation} \tag{ 64 }$$

These states are the archetypal Majorana zero modes [1, 2]. Due to their degeneracy, these modes can be recombined into fermionic quasiparticles by appropriate linear combination, see in equations (4a) and (4b) from section 2.1.

The composition of the eigenvectors slightly changes for the odd case compared to the even N case

$$\begin{equation}{\to {v}}_{\alpha }={\left({x}_{1},\;{y}_{1},\;{x}_{2},\;{y}_{2},\dots ,\;{x}_{\frac{N-1}{2}},\;{y}_{\frac{N-1}{2}},\;{x}_{\frac{N+1}{2}}\right)}^{\text{T}},\end{equation} \tag{ 65 }$$

$$\begin{equation}{\to {v}}_{\beta }={\left({\mathcal{X}}_{1},\;{\mathcal{Y}}_{1},\;{\mathcal{X}}_{2},\;{\mathcal{Y}}_{2},\dots ,\;{\mathcal{X}}_{\frac{N-1}{2}},\;{\mathcal{Y}}_{\frac{N-1}{2}},\;{\mathcal{X}}_{\frac{N+1}{2}}\right)}^{\text{T}}.\end{equation} \tag{ 66 }$$

Although both odd sized chains share the same spectrum, it is possible to find a linear combination of states which belongs to one chain only. The form of the extended states of the odd chains (Δ ≠ ±t and E_± ≠ 0) does not differ much from the one of the even chain and the entries of ${\to {v}}_{\alpha }$ are

$$\begin{equation}\frac{{y}_{l}}{{x}_{1}}=\frac{{E}_{{\pm}}\left({k}_{n}\right)}{a}\;{T}_{l}^{\text{o}}\left({k}_{n}\right),\end{equation} \tag{ 67 }$$

$$\begin{equation}\frac{{x}_{l}}{{x}_{1}}={T}_{l}^{\text{o}}\left({k}_{n}\right)-\frac{b}{a}\;{T}_{l-1}^{\text{o}}\left({k}_{n}\right),\end{equation} \tag{ 68 }$$

where ${T}_{l}^{\text{o}}$ is

$$\begin{equation}{T}_{l}^{\text{o}}\left({k}_{n}\right){:=}\frac{\mathrm{sin}\left(2\;{k}_{n}d\;l\right)}{\mathrm{sin}\left(2\;{k}_{n}d\right)},\end{equation} \tag{ 69 }$$

with k_nd = nπ/(N + 1) (n = 1, ..., N, n ≠ (N + 1)/2). The exchange of a's and b's leads again to the coefficients for the chain β (see appendix C).

The significant difference between even and odd N lies in the realization of the open boundary condition. Solving $\left({\mathcal{H}}_{\alpha }^{\text{odd}}-{E}_{{\pm}}{\mathbb{1}}_{N}\right){\to {v}}_{\alpha }=\to {0}$ yields now

$$\begin{equation}{y}_{0}=0,\quad {y}_{\frac{N+1}{2}}=0,\end{equation} \tag{ 70 }$$

which leads to the momenta k_n.

Zoom In Zoom Out Reset image size

An SSH-like chain with an odd number of sites hosts only a single zero energy mode, but α and β contribute each with one. We find on subchain α for Δ ≠ ±t

$$\begin{equation}{y}_{l}=0,\quad {x}_{l}={\left(\frac{{\Delta}-t}{{\Delta}+t}\right)}^{l-1}\;{x}_{1},\end{equation} \tag{ 71 }$$

and on subchain β

$$\begin{equation}{\mathcal{Y}}_{l}=0,\quad {\mathcal{X}}_{l}={\left(\frac{{\Delta}+t}{{\Delta}-t}\right)}^{l-1}\;{\mathcal{X}}_{1},\end{equation} \tag{ 72 }$$

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 ${\to {\psi }}^{\alpha },\;{\to {\psi }}^{\beta }$, we proceed like for the even N case. The use of the SSH-like basis from equation (22) and the entries of the state ${\to {\psi }}^{\alpha }$ yield now

$$\begin{equation*}{\psi }_{\alpha }^{{\dagger}}=\frac{1}{{v}_{\alpha }}\left[\sum _{j=1}^{\frac{N+1}{2}}\;{x}_{j}\;{\left({\gamma }_{2j-1}^{A}\right)}^{{\dagger}}+\sum _{j=1}^{\frac{N-1}{2}}\;{y}_{j}\;{\left({\gamma }_{2j}^{B}\right)}^{{\dagger}}\right].\end{equation*}$$

Again, the equations (67), (68) and (70) imply a perfect compensation of the A and B sublattice contributions, yielding ${\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}=0$ for E ≠ 0. The zero energy mode, given by its entries in equation (71), leads to ${\left({\psi }_{\alpha }^{{\dagger}}\right)}^{2}=1/2$.

Further, we find that both zero energy modes ${\to {\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{equation}\xi =2d\begin{cases}{\left\vert \mathrm{ln}\left(\frac{t-{\Delta}}{t+{\Delta}}\right)\right\vert }^{-1}\quad \hfill & \vert t\vert {\geqslant}\vert {\Delta}\vert ,\hfill \\ {\left\vert \mathrm{ln}\left(\frac{{\Delta}-t}{{\Delta}+t}\right)\right\vert }^{-1}\quad \hfill & \vert t\vert {\leqslant}\vert {\Delta}\vert ,\hfill \end{cases}\end{equation} \tag{ 73 }$$

and remembering that the atomic site index of x_l is j = 2l − 1 equation (71) yields for |t| ⩾ |Δ|

$$\begin{align}\hfill {x}_{l}& ={x}_{1}\;{\left(-1\right)}^{l-1}\;{\mathrm{e}}^{-2\left(l-1\right)d/\xi },\hfill \\ \hfill & =\;{x}_{1}\;{\left(-1\right)}^{l-1}\;{\mathrm{e}}^{-\left(j-1\right)d/\xi }.\hfill \end{align} \tag{ 74 }$$

For |t| ⩽ |Δ| the x_l coefficients are given by the same equation without the (−1)^l−1 factor. We have moreover q = ±1/ξ, where q is the imaginary momentum yielding E = 0 in equation (41). Thus the localisation of these states is determined only by t and Δ. In the parameter setting of Δ = t we find:

$$\begin{equation}{\to {v}}_{\alpha ,A}={\left(1,\;0,\;0,\dots ,\;0\right)}^{\text{T}},\end{equation} \tag{ 75 }$$

$$\begin{equation}{\to {v}}_{\beta ,B}={\left(0,\;0,\dots ,\;0,\;1\right)}^{\text{T}},\end{equation} \tag{ 76 }$$

while both states exchange their position for Δ = −t

$$\begin{equation}{\to {v}}_{\alpha ,B}={\left(1,\;0,\;0,\dots ,\;0\right)}^{\text{T}},\end{equation} \tag{ 77 }$$

$$\begin{equation}{\to {v}}_{\beta ,A}={\left(0,\;0,\dots ,\;0,\;1\right)}^{\text{T}}.\end{equation} \tag{ 78 }$$

In the last section we have shown that some of the zero energy eigenstates of the BdG Hamiltonian of the finite Kitaev chain are Majorana zero modes (MZM) by exploiting the statistics of the corresponding operators ψ_α,β. We further corroborate this statement now by recalling that an MZM is defined as an eigenstate of the Hamiltonian $\mathcal{H}$ and of the particle hole symmetry $\mathcal{P}$. The latter acts on an eigenstate ${\to {\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 section 2.3, can be constructed in real space too. Of particular interest is their representation in the SSH-like basis. The antiunitary particle–hole symmetry is

$$\begin{equation}\mathcal{P}=\mathcal{K}\;\mathbb{1}.\end{equation} \tag{ 79 }$$

The time reversal and the chiral symmetry depend on N. If N is even we find

$$\begin{equation}{\mathcal{C}}^{\;\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}}=\left[\begin{bmatrix}{cc}\hfill {\mathbb{1}}_{N/2}\otimes {\tau }_{z}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill -{\mathbb{1}}_{N/2}\otimes {\tau }_{z}\hfill \end{bmatrix}\right],\end{equation} \tag{ 80 }$$

$$\begin{equation}{\mathcal{T}}^{\;\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}}=\left[\begin{bmatrix}{cc}\hfill {\mathbb{1}}_{N/2}\otimes {\tau }_{z}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill -{\mathbb{1}}_{N/2}\otimes {\tau }_{z}\hfill \end{bmatrix}\right]\;\mathcal{K}.\end{equation} \tag{ 81 }$$

The expressions for odd N follow by removing the last line and last column in each diagonal block.

The effect of $\mathcal{P}$ from equation (79) can be seen explicitly if one considers $\mathcal{P}\;{\to {\psi }}^{\alpha }$. For N even and Δ ≠ t the elements x_l, y_l of ${\to {\psi }}^{\alpha }$ are given in equations (51) and (52). Here y_l/x₁ is pure imaginary and x_l/x₁ is real. Hence, ${\to {\psi }}^{\alpha }$ is not an eigenstate of $\mathcal{P}$ since the prefactor to y_l is finite, i.e., E_± ≠ 0. We conclude that in a finite Kitaev chain with even number of sites Majorana zero modes emerge only at the Kitaev points for μ = 0, since the states in equations (61)–(64) are eigenstates of $\mathcal{P}$ as well. In the situation of odd N and μ = 0, the eigenstates given by their elements in equations (71), (72) are Majorana zero energy modes for an appropriate choice of x₁. These states can be delocalised over the entire chain, depending on their decay length ξ, while the case of full localisation is only reached at the the Kitaev points, where the MZM turn into the states given by equations (75)–(78).

The last missing situation is to consider a finite chemical potential μ. For this purpose we use the so called chiral basis ${\hat{{\Psi}}}_{c}{:=}{\left({\gamma }_{1}^{A},\;{\gamma }_{2}^{A},\dots ,\;{\gamma }_{N}^{A},\;{\gamma }_{1}^{B},\;{\gamma }_{2}^{B},\dots ,\;{\gamma }_{N}^{B}\right)}^{\text{T}}$. The Kitaev Hamiltonian transforms via ${\hat{H}}_{\text{KC}}=\frac{1}{2}\;{\hat{{\Psi}}}_{c}^{{\dagger}}\;{\mathcal{H}}_{c}\;{\hat{{\Psi}}}_{c}$, into a block off-diagonal matrix

$$\begin{equation}{\mathcal{H}}_{c}=\left[\begin{bmatrix}{cc}\hfill {0}_{N{\times}N}\hfill & \hfill h\hfill \\ \hfill {h}^{{\dagger}}\hfill & \hfill {0}_{N{\times}N}\hfill \end{bmatrix}\right],\end{equation} \tag{ 82 }$$

because there are no ${\gamma }_{j}^{A}\;{\gamma }_{i}^{A}$ (${\gamma }_{j}^{B}\;{\gamma }_{i}^{B}$) contributions in equation (21). The N × N matrix h is tridiagonal

$$\begin{equation}h=\left[\begin{bmatrix}{ccc}\hfill -i\mu \hfill & \hfill a\hfill & \hfill \hfill \\ \hfill -b\hfill & \hfill -i\mu \hfill & \hfill a\hfill \\ \hfill \hfill & \hfill -b\hfill & \hfill -i\mu \hfill & \hfill a\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -b\hfill & \hfill -i\mu \hfill & \hfill a\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -b\hfill & \hfill -i\mu \hfill \end{bmatrix}\right],\end{equation} \tag{ 83 }$$

since the Kitaev Hamiltonian contains only nearest neighbour hoppings. Then the characteristic polynomial is [34]

$$\begin{equation}{P}_{\lambda }\left({\mathcal{H}}_{c}\right)=\mathrm{det}\left({\lambda }^{2}\;{\mathbb{1}}_{N}-h\;{h}^{{\dagger}}\right),\end{equation} \tag{ 84 }$$

where, however, h and h^† do not commute except for t = 0 or Δ = 0. Thus, such matrices cannot be diagonalised simultaneously. Nevertheless the eigenvalues η_j (${\eta }_{j}^{{\ast}}$) of h (h^†) are easily derived e.g. following [33]. We find

$$\begin{equation}{\eta }_{j}=-i\mu +2\;\sqrt{{{\Delta}}^{2}-{t}^{2}}\;\mathrm{cos}\left(\frac{j\pi }{N+1}\right),\quad j=1,\dots ,\;N,\end{equation} \tag{ 85 }$$

independent of whether Δ ⩾ t or t > Δ.

5.1.1. Condition for zero energy modes

Equation (85) immediately yields the criterion for hosting zero energy modes. According to equation (82), we have

$$\begin{equation}\mathrm{det}\left({\mathcal{H}}_{c}\right)=\mathrm{det}\left(h\right)\;\mathrm{det}\left({h}^{{\dagger}}\right)=\vert \mathrm{det}\left(h\right){\vert }^{2},\end{equation} \tag{ 86 }$$

and we need only to focus on det(h).

If a single eigenvalue η_j of h is zero then det(h) vanishes. Thus, for a zero energy mode the chemical potential must satisfy

$$\begin{equation}{\mu }_{n}\;=2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{cos}\left(\frac{n\pi }{N+1}\right),\quad n=1,\dots ,\;N.\end{equation} \tag{ 87 }$$

Obviously, equation (87) cannot be satisfied for generic values of t² − Δ², because all other quantities are real. The only possibility is t² ⩾ Δ². There is only one exception for odd N, because the value n = (N + 1)/2 leads to μ = 0 in equation (87) for all values of t and Δ, in agreement with our results of section 3. This result is exact and confirms findings from [21, 22]; further it improves a similar but approximate condition on the chemical potential discussed by Zvyagin in [23].

An illustration of these discrete solutions μ_n, which we dub 'Majorana lines', is shown in figure 8. All paths contain the Kitaev points at μ = 0 and t = ±Δ. Further, their density is larger close to the boundary of the topological phase, as a result of the slow changes of the cosine function around 0 and π.

Zoom In Zoom Out Reset image size

For growing number of sites N, the density of solutions increases. In the limit N → ∞, θ_n = nπ/N + 1 takes all values in [0, π] and the entire area between $\mu ={\pm}2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}$ for t² ⩾ Δ² 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.

5.1.2. The complete spectrum of the finite Kitaev chain

To proceed we transform equation (84) into an eigenvector problem for hh^†,

$$\begin{equation}h{h}^{{\dagger}}\to {v}={\lambda }^{2}\to {v},\end{equation} \tag{ 88 }$$

where we defined $\to {v}={\left({\xi }_{1},\;{\xi }_{2},\dots ,\;{\xi }_{N}\right)}^{\text{T}}$. Notice that we are not really interested in the eigenvector $\to {v}$ here; we simply use its entries as dummy variables to release a structure hidden in the product of h and h^†. The elements of h,

$$\begin{equation*}{h}_{n,m}=-i\mu \;{\delta }_{n,m}+a\;{\delta }_{n,m+1}-b\;{\delta }_{n+1,m},\end{equation*}$$

where n, m = 1, ..., N, allow us to calculate the product hh^† entry wise

$$\begin{align*}\hfill {\left(h{h}^{{\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]\hfill \\ \hfill & \quad +\;i\mu \;\left(a-b\right)\;\left[{\delta }_{n,m+1}+{\delta }_{n+1,m}\right]\hfill \\ \hfill & \quad +\;ab\;\left({\delta }_{n,m+2}+{\delta }_{n+2,m}\right).\hfill \end{align*}$$

Thus, importantly, equation (88) reveals a recursion formula

$$\begin{align}\hfill {\xi }_{j+2}& =\frac{{\lambda }^{2}+{a}^{2}+{b}^{2}-{\mu }^{2}}{ab}\;{\xi }_{j}-{\xi }_{j-2}-i\mu \left(\frac{a-b}{ab}\right)\hfill \\ \hfill & \hfill \quad {\times}\left({\xi }_{j+1}+{\xi }_{j-1}\right),\end{align} \tag{ 89 }$$

for the components of $\to {v}$. The entries ξ are a generalisation of the Fibonacci polynomials ζ_j from equation (27), to which they reduce for μ = 0, and may be called Tetranacci polynomials [41, 42]. Further, we find the open boundary conditions from equation (88) to be

$$\begin{equation}{\xi }_{0}\;=\;{\xi }_{N+1}=b\;{\xi }_{N+2}-a\;{\xi }_{N}=b\;{\xi }_{1}-a\;{\xi }_{-1}=0,\end{equation} \tag{ 90 }$$

where we used equation (89) for simplifications.

Appendix D contains the description of how to deal with those polynomials, the boundary conditions and further the connection of equation (89) to Kitaev's bulk spectrum λ = E_±(k) in equation (7). Essentially one has to use similar techniques as it was done for the Fibonacci polynomials, where now the power law ansatz ξ_j ∝ r^j leads to a characteristic equation for r of order four. Thus, we find in total four linearly independent fundamental solutions r_±1,±2, which can be expressed in terms of two complex wavevectors denoted by κ_1,2 through the equality

$$\begin{equation}{r}_{{\pm}j}={\mathrm{e}}^{{\pm}\text{i}{\kappa }_{j}},\quad j=1,2.\end{equation} \tag{ 91 }$$

These wavevectors are not independent, but coupled via

$$\begin{equation}\mathrm{cos}\left({\kappa }_{1}\right)+\mathrm{cos}\left({\kappa }_{2}\right)=-\frac{\mu t}{{t}^{2}-{{\Delta}}^{2}},\quad \forall \;t,\;{\Delta},\;\mu \in \mathbb{R}.\end{equation} \tag{ 92 }$$

For μ = 0 we can recover from equation (92) our previous results, whereby one has only pure real (k) or pure imaginary (iq) wavevectors⁵. Further, equations (7) and (92) yield

$$\begin{equation*}{E}_{{\pm}}\left({\kappa }_{1}\right)={E}_{{\pm}}\left({\kappa }_{2}\right).\end{equation*}$$

The linearity of the recursion formula equation (89) states that the superposition of all four fundamental solutions is the general form of ξ_j. Since the boundary conditions translate into a homogeneous system of four coupled equations and a trivial solution for ξ_j has to be avoided, we find that the determinant of the matrix describing these equations has to be zero. After some algebraic manipulations, this procedure leads finally to the full quantization rule of the Kitaev chain

$$\begin{equation}F\left({\kappa }_{1},{\kappa }_{2}\right)=F\left({\kappa }_{1},-{\kappa }_{2}\right),\end{equation} \tag{ 93 }$$

where we introduced the function F(κ₁, κ₂) as

$$\begin{align}\hfill F\left({\kappa }_{1},{\kappa }_{2}\right)& ={\mathrm{sin}}^{2}\left[\frac{{\kappa }_{1}+{\kappa }_{2}}{2}\;\left(N+1\right)\right]\hfill \\ \hfill & \hfill \quad {\times}\left[1+{\left(\frac{{\Delta}}{t}\right)}^{2}\;{\mathrm{cot}}^{2}\left(\frac{{\kappa }_{1}+{\kappa }_{2}}{2}\right)\right].\end{align} \tag{ 94 }$$

Similar quantization conditions are known for an open X–Y spin chain in transverse field [26]. Notice that the quantization rule is symmetric with respect to κ_1,2. Table 1 gives an overview of the quantization rules for different parameter settings $\left({\Delta},\;t,\;\mu \right)$. The bulk eigenvalues of a finite Kitaev chain with four sites and μ ≠ 0 are shown in figure 9.

Zoom In Zoom Out Reset image size

The previous relations open another route to finding the condition leading to modes with exact zero energy. A convenient form of equation (92) is

$$\begin{equation}\mathrm{cos}\left(\frac{{\kappa }_{1}+{\kappa }_{2}}{2}\right)\;\mathrm{cos}\left(\frac{{\kappa }_{1}-{\kappa }_{2}}{2}\right)=-\frac{1}{2}\;\frac{\mu t}{{t}^{2}-{{\Delta}}^{2}},\end{equation} \tag{ 95 }$$

and the dispersion relation can be transformed into

$$\begin{align}\hfill {E}^{2}& =\frac{1}{{\mathrm{cos}}^{2}\left(\frac{{\kappa }_{1}{\pm}{\kappa }_{2}}{2}\right)}\left[4\left({t}^{2}-{{\Delta}}^{2}\right){\mathrm{cos}}^{2}\left(\frac{{\kappa }_{1}{\pm}{\kappa }_{2}}{2}\right)-{\mu }^{2}\right]\hfill \\ \hfill & \quad {\times}\;\left[\frac{{t}^{2}}{{t}^{2}-{{\Delta}}^{2}}-{\mathrm{cos}}^{2}\left(\frac{{\kappa }_{1}{\pm}{\kappa }_{2}}{2}\right)\right].\hfill \end{align} \tag{ 96 }$$

Both combinations κ₁ ± κ₂ yield the same energy, due to equation (95). If one of the brackets in equation (96) vanishes for κ₁ + κ₂ (κ₁ − κ₂), the second one does so for κ₁ − κ₂ (κ₁ + κ₂) too. Hence, zero energy is achieved exactly if

$$\begin{equation}1+{\left(\frac{{\Delta}}{t}\right)}^{2}\;{\mathrm{cot}}^{2}\left(\frac{{\kappa }_{1}{\pm}{\kappa }_{2}}{2}\right)=0.\end{equation} \tag{ 97 }$$

This puts restrictions on Δ, t, μ. Together with the quantization rule in equation (93), this ultimately leads to (87) 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 → ∞ the difference between even and odd N vanishes and the part of the μ = 0 axis between t = −Δ and t = Δ for even N leads to zero energy states too, in virtue of equation (42). Analogously, the area around the origin in figure 8, defined by $2\sqrt{{t}^{2}-{{\Delta}}^{2}}{< }\vert \mu \vert$ and μ < 2|t| with μ ≠ 0 does not support zero energy modes for all finite N. Instead, this area contains solutions with exponentially small energies, see equation (96), which become zero exclusively in the limit N → ∞. The wavevectors obey κ₁ = iq₁, κ₂ = π + iq₂ with real q_1,2 for Δ² > t², and κ₁ = iq₁, κ₂ = iq₂ otherwise, which follows from equation (93) after some manipulations; see also [26]. Thus, the entire non trivial phase hosts zero energy solutions for N → ∞.

The calculation of an arbitrary wave function of the Kitaev Hamiltonian without any restriction on t, Δ, μ is performed at best in the chiral basis yielding the block off-diagonal structure in equation (82). A suitable starting point is to consider a vector $\to {w}$ in the following form

$$\begin{equation*}\to {w}=\left(\begin{pmatrix}\hfill \to {v}\hfill \\ \hfill \to {u}\hfill \end{pmatrix}\right)\end{equation*}$$

with $\to {v}=\left({\xi }_{1},\dots ,\;{\xi }_{N}\right)$, $\to {u}=\left({\sigma }_{1},\dots ,\;{\sigma }_{N}\right)$. Solving for an eigenstate with eigenvalue λ demands on $\to {v},\;\to {u}$

$$\begin{equation}h\;\to {u}=\lambda \;\to {v},\end{equation} \tag{ 98 }$$

$$\begin{equation}{h}^{{\dagger}}\;\to {v}=\lambda \;\to {u},\end{equation} \tag{ 99 }$$

with h from equation (83). Thus, $h\;{h}^{{\dagger}}\to {v}={\lambda }^{2}\;\to {v}$, and we recover equation (88) and the entries of $\to {v}$ obey equation (89) again. In appendix E we derived the closed formula for ξ_j, namely

$$\begin{equation}{\xi }_{j}=\sum _{i=-2}^{1}\;{\xi }_{i}\;{X}_{i}\left(j\right),\quad j\in \mathbb{Z}\end{equation} \tag{ 100 }$$

where ξ₋₂, ..., ξ₁ are the initial values of the polynomial sequence dependent on the boundary conditions, and X_i(j) inherit the selective property

$$\begin{equation}{X}_{i}\left(j\right)={\delta }_{i,j}\quad \mathrm{f}\mathrm{o}\mathrm{r}\;\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{y}i,\;j=-2,\dots ,\;1.\end{equation} \tag{ 101 }$$

That these functions X_i(j) exist and that they indeed satisfy equation (100) for arbitrary values of j is discussed in appendix E.

The remaining task is to obtain the initial values ξ₋₂, ..., ξ₁, which follow from the open boundary conditions equation (90). Further, one has one free degree of freedom, which we to choose to be the entry ξ₁ of $\to {v}$. In total our initial values are ξ₁, ξ₀ = 0; ξ₋₁ = bξ₁/a and ξ₋₂ follows from ξ_N+1 = 0 and equation (100),

$$\begin{equation*}{\xi }_{-2}=-{\xi }_{1}\;\frac{a\;{X}_{1}\left(N+1\right)+b\;{X}_{-1}\left(N+1\right)}{a\;{X}_{-2}\left(N+1\right)}.\end{equation*}$$

Demanding further bξ_N+2 − aξ_N = 0 quantizes the momenta κ_1,2 and in turn λ = E_±(κ_1,2), according to equation (93) (the relation between X_j and κ_1,2 is discussed in the appendix E). Notice that the form given by equations (100) and (102) (below) hold for all eigenstates of the Kitaev BdG Hamiltonian; the distinction between extended/decaying states and MZM is made by the values of k_1,2, or equivalently κ_1,2.

The second part of the eigenstate $\to {w}$, i.e. $\to {u}$, follows in principle from equation (99). A simpler and faster way is to consider ${h}^{{\dagger}}\;h\;\to {u}={\lambda }^{2}\;\to {u}$. A comparison of h and h^† reveals that they transform into each other by exchanging a and b and switching μ into −μ. Consequently, the structure of the entries of $\to {u}$ follow essentially from the ones of $\to {v}$. Thus, we find that the σ_i obey equation (89) too, yielding

$$\begin{equation}{\sigma }_{j}=\sum _{i=-2}^{1}\;{\sigma }_{i}\;{X}_{i}\left(j\right),\quad j\in \mathbb{Z}\end{equation} \tag{ 102 }$$

with the same functions X_i(j). The boundary condition on $\to {u}$ reads

$$\begin{equation*}{\sigma }_{0}\;=\;{\sigma }_{N+1}=a\;{\sigma }_{N+2}-b\;{\sigma }_{N}=a\;{\sigma }_{1}-b\;{\sigma }_{-1}=0.\end{equation*}$$

Proceeding as we did for $\to {v}$ yields σ₀ = 0, σ₋₁ = aσ₁/b and

$$\begin{equation}{\sigma }_{-2}=-{\sigma }_{1}\;\frac{b\;{X}_{1}\left(N+1\right)+a\;{X}_{-1}\left(N+1\right)}{b\;{X}_{-2}\left(N+1\right)},\end{equation} \tag{ 103 }$$

where σ₁ is fixed by the first line of equation (99)

$$\begin{equation}{\sigma }_{1}=\frac{i\mu \;{\xi }_{1}+b\;{\xi }_{2}}{\lambda }.\end{equation} \tag{ 104 }$$

The last open boundary condition aσ_N+2 − bσ_N = 0 is satisfied for λ = E_±(κ_1,2) and thus already by the construction of $\to {v}$. We notice that we assumed λ ≠ 0 in order to obtain σ₁; the case λ = 0 is discussed in section 4 below.

In the limit of μ = 0 on $\to {w}$ it holds that

$$\begin{align*}\hfill {X}_{-2}\left(2l+1\right){\vert }_{\mu =0}& =0,\hfill \\ \hfill {X}_{0}\left(2l+1\right){\vert }_{\mu =0}& =0,\hfill \\ \hfill {X}_{-1}\left(2l\right){\vert }_{\mu =0}& =0,\hfill \\ \hfill {X}_{1}\left(2l\right){\vert }_{\mu =0}& =0,\hfill \end{align*}$$

for all values of l; thus only two initial values for ξ_j and two for σ_j are necessary to fix the sequences.

Finally, one can prove easily that the functions X_i(j) are always real and in consequence all ξ_j (σ_j) are real (pure imaginary) if ξ₁ is chosen to be real. Thus, the corresponding operators ${\psi }_{\to {w}}\;\left({\psi }_{\to {w}}^{{\dagger}}\right)$ can never square to 1/2.

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The full analytic calculation of the spectrum is at best performed in the basis of Majorana operators ${\gamma }_{j}^{A\left(B\right)}$, ordered according to the chain index

$$\begin{equation}{\hat{{\Psi}}}_{M,\;co}\;{:=}\;{\left({\gamma }_{1}^{A},\;{\gamma }_{1}^{B},\;{\gamma }_{2}^{A},\;{\gamma }_{2}^{B},\dots ,\;{\gamma }_{N}^{A},\;{\gamma }_{N}^{B}\right)}^{\text{T}}.\end{equation} \tag{ B1 }$$

Then the BdG Hamiltonian becomes block tridiagonal

$$\begin{equation}{\mathcal{H}}_{M,\text{co}}^{\text{KC}}={\left[\begin{bmatrix}{cc}\hfill A\hfill & \hfill B\hfill \\ \hfill {B}^{{\dagger}}\hfill & \hfill A\hfill & \hfill B\hfill \\ \hfill \hfill & \hfill {B}^{{\dagger}}\hfill & \hfill A\hfill & \hfill B\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill {B}^{{\dagger}}\hfill & \hfill A\hfill & \hfill B\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill {B}^{{\dagger}}\hfill & \hfill A\hfill \end{bmatrix}\right]}_{2N{\times}2N},\end{equation} \tag{ B2 }$$

where A and B are 2 × 2 matrices

$$\begin{equation}A=\left[\begin{bmatrix}{cc}\hfill 0\hfill & \hfill -i\mu \hfill \\ \hfill i\mu \hfill & \hfill 0\hfill \end{bmatrix}\right],\quad B=\left[\begin{bmatrix}{cc}\hfill 0\hfill & \hfill a\hfill \\ \hfill b\hfill & \hfill 0\hfill \end{bmatrix}\right].\end{equation} \tag{ B3 }$$

Since we are interested in the spectrum, we have essentially only to calculate (and factorise) the characteristic polynomial ${P}_{\lambda }\left({\mathcal{H}}^{\text{KC}}\right)=\mathrm{det}\left(\lambda \mathbb{1}-{\mathcal{H}}_{\text{KC}}\right)$ which reads simply

$$\begin{equation}{P}_{\lambda }=\mathrm{det}{\left[\begin{bmatrix}{cc}\hfill \lambda \mathbb{1}\hfill & \hfill -B\hfill \\ \hfill -{B}^{{\dagger}}\hfill & \hfill \lambda \mathbb{1}\hfill & \hfill -B\hfill \\ \hfill \hfill & \hfill -{B}^{{\dagger}}\hfill & \hfill \lambda \mathbb{1}\hfill & \hfill -B\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -{B}^{{\dagger}}\hfill & \hfill \lambda \mathbb{1}\hfill & \hfill -B\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -{B}^{{\dagger}}\hfill & \hfill \lambda \mathbb{1}\hfill \end{bmatrix}\right]}_{2N{\times}2N},\end{equation} \tag{ B4 }$$

at zero μ. In the following we will consider λ to be just a 'parameter', which is not necessarily real in the beginning. Further, we shall impose (only in the beginning) the restriction λ ≠ 0. However, our results will even hold without them. The validity of our argument follows from the fact that the determinant P_λ is a smooth function

$$\begin{equation}{P}_{\lambda }\in {C}^{\infty }\left(\mathbb{P}\right),\end{equation} \tag{ B5 }$$

in the entire parameter space $\mathbb{P}{:=}{\mathbb{C}}^{3}$, which contains a, b and λ.

The technique we want to use to evaluate P_λ is essentially given by the recursion formula of the 2 × 2 matrices Λ_j [43, 44]:

$$\begin{equation}{{\Lambda}}_{j}=\lambda {\mathbb{1}}_{2}-{B}^{{\dagger}}{{\Lambda}}_{j-1}^{-1}B,\quad {{\Lambda}}_{1}{:=}\lambda {\mathbb{1}}_{2},\end{equation} \tag{ B6 }$$

where j = 1, ..., N and ${P}_{\lambda }=\prod _{j=1}^{N}\mathrm{det}\left({{\Lambda}}_{j}\right)$.

The matrices B and B^† are pure off-diagonal matrices and since $\lambda {\mathbb{1}}_{2}$ is diagonal, one can prove that Λ_j has the general diagonal form of ${{\Lambda}}_{j}{:=}\left[\begin{bmatrix}{cc}\hfill {x}_{j}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {y}_{j}\hfill \end{bmatrix}\right]$ (for all j). The application of equation (B6) leads to a recursion formula for both sequences of entries

$$\begin{align*}\hfill {x}_{j+1}& =\lambda +\frac{{b}^{2}}{{y}_{j}},\hfill \\ \hfill {y}_{j+1}& =\lambda +\frac{{a}^{2}}{{x}_{j}},\hfill \end{align*}$$

and the initial values are x₁ = y₁ = λ. We find x_j and y_j to be fractions in general, and define ζ_j, _j, β_j and δ_j by

$$\begin{align*}\hfill {x}_{j}& =:\;\frac{{\zeta }_{j}}{{\beta }_{j}},\hfill \\ \hfill {y}_{j}& =:\;\frac{{{\epsilon}}_{j}}{{\delta }_{j}},\hfill \end{align*}$$

to take this into account. The initial values can be set as

$$\begin{equation}{\zeta }_{1}={{\epsilon}}_{1}=\lambda ,\end{equation} \tag{ B7 }$$

$$\begin{equation}{\beta }_{1}={\delta }_{1}=1,\end{equation} \tag{ B8 }$$

and after a little bit of algebra we find their growing rules to be

$$\begin{equation}{\zeta }_{j+1}=\lambda \;{{\epsilon}}_{j}+{b}^{2}\;{\delta }_{j},\end{equation} \tag{ B9 }$$

$$\begin{equation}{{\epsilon}}_{j+1}=\lambda \;{\zeta }_{j}+{a}^{2}\;{\beta }_{j},\end{equation} \tag{ B10 }$$

$$\begin{equation}{\beta }_{j+1}={{\epsilon}}_{j},\end{equation} \tag{ B11 }$$

$$\begin{equation}{\delta }_{j+1}={\zeta }_{j},\end{equation} \tag{ B12 }$$

where j starts from 1. The definitions ζ₀ := δ₁ = 1 and ₀ := β₁ = 1, enable us to get rid of the δ_j and β_j terms inside equations (B9) and (B10). Hence

$$\begin{equation}{\zeta }_{j+1}=\lambda \;{{\epsilon}}_{j}+{b}^{2}\;{\zeta }_{j-1},\end{equation} \tag{ B13 }$$

$$\begin{equation}{{\epsilon}}_{j+1}=\lambda \;{\zeta }_{j}+{a}^{2}\;{{\epsilon}}_{j-1},\end{equation} \tag{ B14 }$$

which leads to the relations

$$\begin{equation}{\zeta }_{2}={\lambda }^{2}+{b}^{2},\end{equation} \tag{ B15 }$$

$$\begin{equation}{{\epsilon}}_{2}={\lambda }^{2}+{a}^{2}.\end{equation} \tag{ B16 }$$

We already extended the sequences of ζ_j and _j artificially backwards and we continue to do so, using the equations (B13) and (B14), starting from j = −1 with ζ₋₁ = ₋₁ = 0. Please note there are no corresponding x₀, y₀ or even x₋₁, y₋₁ expressions, since they would involve division by 0.

The last duty of β_j and δ_j is to simplify the determinant P_λ by using the equations (B8), (B11) and (B12)

$$\begin{equation}{P}_{\lambda }=\prod _{j=1}^{N}\mathrm{det}\left({{\Lambda}}_{j}\right)=\prod _{j=1}^{N}{x}_{j}\;{y}_{j}={\zeta }_{N}\;{{\epsilon}}_{N},\end{equation} \tag{ B17 }$$

which reduces the problem to finding only ζ_N and _N.

Please note that the determinant is in fact independent of the choice of the initial values for ζ₁, ₁, β₁ and δ₁ in the equations (B7) and (B8). Further, equations (B13), (B14) and (B17) together show the predicted smoothness of P_λ in $\mathbb{P}$ and all earlier restrictions are not important anymore. Finally we consider λ to be real again.

Even though it seems that we are left with the calculation of two polynomials, we need in fact only one, because both are linked via the exchange of a and b. Note that λ is considered here as a number and thus does not depend on a and b. Further, the dispersion relation is invariant under this exchange.

The connection of ζ_j and _j for all j ⩾ −1 is

$$\begin{align*}\hfill {\zeta }_{j}\equiv {\zeta }_{j}\left(a,b\right)& ={{\epsilon}}_{j}\left(b,a\right),\hfill \\ \hfill {{\epsilon}}_{j}\equiv {{\epsilon}}_{j}\left(a,b\right)& ={\zeta }_{j}\left(b,a\right),\hfill \end{align*}$$

and can be proven via induction using equations (B13) and (B14). Decoupling ζ_j and _j yields

$$\begin{equation}{\zeta }_{j+2}=\left[{\lambda }^{2}+{a}^{2}+{b}^{2}\right]\;{\zeta }_{j}-{a}^{2}{b}^{2}\;{\zeta }_{j-2},\end{equation} \tag{ B18 }$$

where one identifies them as (generalized) Fibonacci polynomials [36, 37]. The qualitative difference between even and odd number of sites is a consequence of equation (B18) and the initial values for ζ_j.

The next step is to obtain the closed form expression of ζ_j (_j), the so called Binet form. We focus exclusively on ζ_j.

One way to keep the notation easier is to introduce x := λ² + a² + b², y := a²b², v_j := ζ_2j and u_j := ζ_2j−1, such that u_j (v_j) obey

$$\begin{equation*}{u}_{j+1}=x\;{u}_{j}-y\;{u}_{j-1}.\end{equation*}$$

The Binet form can be obtained by using a power law ansatz u_j ∝ r^j, leading to two fundamental solutions

$$\begin{equation}{r}_{1,2}=\frac{x\;{\pm}\;\sqrt{{x}^{2}-4\;y}}{2}.\end{equation} \tag{ B19 }$$

Please note that this square root is always well defined, which can be seen in the simplest way by setting λ to zero. Consequently, the difference between r₁ and r₂ is never zero.

A general solution of u_j (v_j) can be achieved with a superposition of r_1,2 with some coefficients c_1,2,

$$\begin{equation*}{u}_{j}={c}_{1}\;{r}_{1}^{j}+{c}_{2}\;{r}_{2}^{j},\end{equation*}$$

due to the linearity of their recursion formula. Both constants c₁ and c₂ are fixed by the initial values for ζ_j, for example u₀ = ζ₋₁ = 0 and u₁ = ζ₁ = λ and similar for v_j. After some simplifications, we finally arrive at

$$\begin{equation}{\zeta }_{2j-1}=\lambda \;\frac{{r}_{1}^{j}-{r}_{2}^{j}}{{r}_{1}-{r}_{2}},\end{equation} \tag{ B20 }$$

$$\begin{equation}{\zeta }_{2j}=\frac{\left[{\lambda }^{2}+{b}^{2}\right]\;\left({r}_{1}^{j}-{r}_{2}^{j}\right)-{r}_{1}\;{r}_{2}\left({r}_{1}^{j-1}-{r}_{2}^{j-1}\right)}{{r}_{1}-{r}_{2}},\end{equation} \tag{ B21 }$$

in agreement with [36–38]. The validity of the solutions is guaranteed by a proof via induction, where one needs mostly the properties of r_1,2 to be the fundamental solutions. The exchange of a and b leads to the expressions

$$\begin{equation}{{\epsilon}}_{2j-1}=\lambda \;\frac{{r}_{1}^{j}-{r}_{2}^{j}}{{r}_{1}-{r}_{2}},\end{equation} \tag{ B22 }$$

$$\begin{equation}{{\epsilon}}_{2j}=\frac{\left[{\lambda }^{2}+{a}^{2}\right]\;\left({r}_{1}^{j}-{r}_{2}^{j}\right)-{r}_{1}\;{r}_{2}\left({r}_{1}^{j-1}-{r}_{2}^{j-1}\right)}{{r}_{1}-{r}_{2}},\end{equation} \tag{ B23 }$$

where we used that r_1,2 is symmetric in a and b. At this stage we have the characteristic polynomial in closed form for all Δ, t and more importantly for all sizes N at zero μ.

We can already anticipate the twice degenerated eigenvalues of the odd sized Kitaev chain, because from the closed forms of _j and ζ_j it follows immediately

$$\begin{equation}{\zeta }_{\text{odd}}={{\epsilon}}_{\text{odd}}.\end{equation} \tag{ B24 }$$

Notice that equation (B24) is important to derive the characteristic polynomial via the SSH description of the Kitaev BdG Hamiltonian at μ = 0 and to show the equivalence to the approach used here. It is recommended to use the determinant formula in [45] together with equations (B9) and (B10) for the proof.

The main steps of the factorisation are mentioned in the next section.

The trick to factorise our Fibonacci polynomials [36, 37] bases on the special form of r_1,2. The ansatz is to look for the eigenvalues λ in the following form

$$\begin{equation}x=2\sqrt{y}\;\mathrm{cos}\left(\theta \right),\end{equation} \tag{ B25 }$$

which is actually the definition of θ. The hermiticity of the Hamiltonian enforces real eigenvalues and consequently θ can be chosen either real, describing extended solutions, or pure imaginary, which is connected to decaying states. The ansatz leads to an exponential form of the fundamental solutions

$$\begin{align*}\hfill {r}_{1}& =\sqrt{y}\;{\mathrm{e}}^{\text{i}\;\theta },\hfill \\ \hfill {r}_{2}& =\sqrt{y}\;{\mathrm{e}}^{-\text{i}\;\theta },\hfill \end{align*}$$

and we consider $\theta \in \mathbb{R}$ first. Thus, we find the eigenvalues for odd N

$$\begin{equation*}{{\epsilon}}_{N}={\zeta }_{N}=\lambda \;\frac{\mathrm{sin}\left(\frac{N+1}{2}\;\theta \right)}{\mathrm{sin}\left(\theta \right)}\;{\sqrt{y}}^{\frac{N-1}{2}}\;=\;0.\end{equation*}$$

One obvious solution is λ = 0. The introduction of 2kd = θ, where d is the lattice constant of the Kitaev chain, leads to:

$$\begin{equation}\frac{\mathrm{sin}\left[\left(N+1\right)\;kd\right]}{\mathrm{sin}\left(2\;kd\right)}=\;0,\end{equation} \tag{ B26 }$$

and solutions inside the first Brillouin-zone are given by

$$\begin{equation*}{k}_{n}d=\frac{n\pi }{N+1}\end{equation*}$$

where n runs from 1, ..., N without (N + 1)/2. Please note that equation (B26) cannot be satisfied for N = 1.

The even N case requires more manipulations. We first rearrange equation (B21) as

$$\begin{equation*}{\zeta }_{2j}=\frac{\left({\lambda }^{2}+{b}^{2}-{r}_{2}\right)\;{r}_{1}^{j}-\left({\lambda }^{2}+{b}^{2}-{r}_{1}\right)\;{r}_{2}^{j}}{{r}_{1}-{r}_{2}}.\end{equation*}$$

The expressions λ² + b² − r_1,2 are simplified to

$$\begin{align*}\hfill {\lambda }^{2}+{b}^{2}-{r}_{1}& =x-{a}^{2}-{r}_{1}=\sqrt{y}\;{\mathrm{e}}^{-\text{i}\theta }-{a}^{2},\hfill \\ \hfill {\lambda }^{2}+{b}^{2}-{r}_{2}& =\sqrt{y}\;{\mathrm{e}}^{\text{i}\theta }-{a}^{2}.\hfill \end{align*}$$

In the end ζ_2j becomes

$$\begin{equation}{\zeta }_{2j}={\left(\sqrt{y}\right)}^{j-1}\frac{\sqrt{y}\;\mathrm{sin}\left[\theta \left(j+1\right)\right]-{a}^{2}\;\mathrm{sin}\left(\theta j\right)}{\mathrm{sin}\left(\theta \right)}.\end{equation} \tag{ B27 }$$

Note that the competition of Δ and t is hidden inside the square root

$$\begin{equation*}\sqrt{y}=\left\{\begin{matrix}\hfill {{\Delta}}^{2}-{t}^{2},\hfill & \hfill \mathrm{i}\mathrm{f}\;\vert {\Delta}\vert { >}\vert t\vert \hfill \\ \hfill {t}^{2}-{{\Delta}}^{2},\hfill & \hfill \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}\hfill \end{matrix}\right.,\end{equation*}$$

affecting both the quantization condition and the dispersion relation E_±(k) = λ(θ), which follow from equation (B25). However, both situations lead to the same result, because the momenta and the spectrum are shifted by π/2 (with respect to kd). From ζ_N it follows:

$$\begin{equation*}{\Delta}\;\frac{\mathrm{sin}\left[kd\;\left(N+1\right)\right]\;\mathrm{cos}\left(kd\right)}{\mathrm{sin}\left(2\;kd\right)}+t\;\frac{\mathrm{cos}\left[kd\;\left(N+1\right)\right]\;\mathrm{sin}\left(kd\right)}{\mathrm{sin}\left(2\;kd\right)}=0,\end{equation*}$$

or in shorter form

$$\begin{equation}\mathrm{tan}\left[kd\;\left(N+1\right)\right]=-\frac{{\Delta}}{t}\mathrm{tan}\left(kd\right),\quad kd\ne 0,\frac{\pi }{2}\end{equation} \tag{ B28 }$$

for even N. The polynomial _N can be treated in the same way leading to

$$\begin{equation}\mathrm{tan}\left[kd\;\left(N+1\right)\right]=\frac{{\Delta}}{t}\mathrm{tan}\left(kd\right),\quad kd\ne 0,\frac{\pi }{2}.\end{equation} \tag{ B29 }$$

From equation (B25) follows the bulk spectrum for all N,

$$\begin{equation*}\lambda \left(\theta \right)={E}_{{\pm}}\left(kd\right)={\pm}\sqrt{4{{\Delta}}^{2}\;{\mathrm{sin}}^{2}\left(kd\right)+4{t}^{2}\;{\mathrm{cos}}^{2}\left(kd\right)},\end{equation*}$$

in agreement with equations (34) and (35).

The case of decaying states is similar, but not just done by replacing k by iq. The following case is only valid for even N, since we have already all 2N eigenvalues of the odd N case.

Our ansatz is modified to

$$\begin{equation*}x=-2\sqrt{y}\;\mathrm{cosh}\left(\theta \right),\end{equation*}$$

by an additional minus sign, which is important to find the decaying state solutions. After some manipulations ζ_N = 0 yields the quantization conditions

$$\begin{equation}\frac{t\;\mathrm{tanh}\left[qd\;\left(N+1\right)\right]+{\Delta}\;\mathrm{tanh}\left(qd\right)}{\mathrm{sinh}\left(2\;qd\right)}=0,\quad \vert {\Delta}\vert {\geqslant}\vert t\vert ,\end{equation} \tag{ B30 }$$

$$\begin{equation}\frac{{\Delta}\;\mathrm{tanh}\left[qd\;\left(N+1\right)\right]+t\;\mathrm{tanh}\left(qd\right)}{\mathrm{sinh}\left(2\;qd\right)}=0,\quad \vert t\vert {\geqslant}\vert {\Delta}\vert ,\end{equation} \tag{ B31 }$$

where qd = θ/2. The conditions for qd = 0 as solution, corresponding to infinite decay length ξ, turn out to be ±t/Δ = N + 1 (if |t| ⩾ |Δ|) or ±Δ/t = N + 1 (else) and follow by applying the limit qd → 0 on equations (B30) and (B31).

A last simplification can be done for qd ≠ 0

$$\begin{equation*}\mathrm{tanh}\left[qd\;\left(N+1\right)\right]=-m\;\mathrm{tanh}\left(qd\right),\end{equation*}$$

where we introduced

$$\begin{equation*}m=\left\{\begin{matrix}\hfill \frac{{\Delta}}{t}\hfill & \hfill \mathrm{i}\mathrm{f}\vert {\Delta}\vert { >}\vert t\vert \hfill \\ \hfill \frac{t}{{\Delta}}\hfill & \hfill \mathrm{i}\mathrm{f}\vert t\vert { >}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation*}$$

The criterion to find a wave vector is that (−m) ⩾ 1, but not larger than N + 1, which leads then to exactly two solutions ±q and otherwise to none. The corresponding eigenvalues can be obtained from

$$\begin{align*}\hfill {E}_{{\pm}}\left(qd\right)={\pm}\sqrt{4{t}^{2}\;{\mathrm{cosh}}^{2}\left(qd\right)-4{{\Delta}}^{2}\;{\mathrm{sinh}}^{2}\left(qd\right)},\quad \vert {\Delta}\vert {\geqslant}\vert t\vert ,\\ \hfill {E}_{{\pm}}\left(qd\right)={\pm}\sqrt{4{{\Delta}}^{2}\;{\mathrm{cosh}}^{2}\left(qd\right)-4{t}^{2}\;{\mathrm{sinh}}^{2}\left(qd\right)},\quad \vert t\vert {\geqslant}\vert {\Delta}\vert ,\end{align*}$$

which can be zero. The results for _N can be obtained by replacing t with −t everywhere.

The sublattice vectors are defined via the N entries

$$\begin{align*}\hfill {\to {v}}_{\alpha }& ={\left({x}_{1},\;{y}_{1},\;{x}_{2},\;{y}_{2},\dots ,\;{x}_{N/2},\;{y}_{N/2}\right)}^{\text{T}},\hfill \\ \hfill {\to {v}}_{\beta }& ={\left({\mathcal{X}}_{1},\;{\mathcal{Y}}_{1},\;{\mathcal{X}}_{2},\;{\mathcal{Y}}_{2},\dots ,\;{\mathcal{X}}_{N/2},\;{\mathcal{Y}}_{N/2}\right)}^{\text{T}}.\hfill \end{align*}$$

Solving $\left({\mathcal{H}}_{\alpha }^{\text{even}}-\lambda \;\mathbb{1}\right)\;{\to {v}}_{\alpha }=0$ leads to

$$\begin{equation}\;a\;{y}_{1}=\lambda \;{x}_{1},\end{equation} \tag{ C1 }$$

$$\begin{equation}-a\;{x}_{N/2}=\lambda \;{y}_{N/2},\end{equation} \tag{ C2 }$$

and

$$\begin{equation}b\;{x}_{l+1}-a\;{x}_{l}=\lambda \;{y}_{l},\end{equation} \tag{ C3 }$$

$$\begin{equation}a\;{y}_{l+1}-b\;{y}_{l}=\lambda \;{x}_{l+1},\end{equation} \tag{ C4 }$$

where l runs from 1 to (N/2) − 1. The coupled equations (C1)–(C4) for the entries of the eigenvector are continuous in all parameters. However, resolving to the x_l's and y_l's may lead to problems for certain values of Δ, t and λ.

Case 1. |Δ| ≠ |t|. The parameter setting excludes λ = 0, as we found from our spectral analysis in section 3. Equations (C3) and (C4) are used to disentangle x's and y's. Both sequences obey

$$\begin{equation}{y}_{l+1}=\frac{{\lambda }^{2}+{a}^{2}+{b}^{2}}{ab}\;{y}_{l}-{y}_{l-1},\end{equation} \tag{ C5 }$$

where l = 1, ..., N/2. Thus the y's and x's are Fibonacci polynomials [36, 37]. The difference to the previous ones found for the spectrum is that the new version can be dimensionless in physical units, depending on the initial values. The transformation formula [37] to pass from the unitless recursion formula to the other one is given by equation (A6).

The Binet form of the dimensionless sequences is obtained with same treatment as for the spectrum. The power ansatz y_l ∝ f^l yields the fundamental solutions f_1,2, obeying

$$\begin{equation}{f}_{1}+{f}_{2}=\frac{{\lambda }^{2}+{a}^{2}+{b}^{2}}{ab},\end{equation} \tag{ C6 }$$

$$\begin{equation}{f}_{1}\;\cdot \;{f}_{2}=1,\end{equation} \tag{ C7 }$$

$$\begin{equation}{f}_{1}\;\ne \;{f}_{2}.\end{equation} \tag{ C8 }$$

Due to the linearity of the recursion formula, the most generic ansatz for y_l is

$$\begin{equation*}{y}_{l}={c}_{1}\;{f}_{1}^{l}+{c}_{2}\;{f}_{2}^{l},\end{equation*}$$

where the constants c_{1, 2} follow from the initial values of y_{1, 2}. The calculation of both constants leads to

$$\begin{equation*}{y}_{l}={y}_{2}\;{T}_{l-1}-{y}_{1}\;{T}_{l-2},\end{equation*}$$

where T_l is simply [40]

$$\begin{equation*}{T}_{l}{:=}\frac{{f}_{1}^{l}-{f}_{2}^{l}}{{f}_{1}-{f}_{2}}.\end{equation*}$$

Analogously we find

$$\begin{equation*}{x}_{l}={x}_{2}\;{T}_{l-1}-{x}_{1}\;{T}_{l-2}.\end{equation*}$$

A short comment on the initial values y_1,2. A hermitian matrix is always diagonalisable, regardless of degeneracies in its spectrum and an eigenvector is well defined only up to the prefactor. Consequently we have the freedom to choose one component of ${\to {v}}_{\alpha }$. This choice will in turn define all remaining initial values.

Consider for example x₁ to be a fixed value of our choice. We find y₁, x₂ and y₂ to be

$$\begin{align*}\hfill {y}_{1}& =\frac{\lambda }{a}\;{x}_{1},\hfill \\ \hfill {y}_{2}& =\frac{\lambda }{a}\;\frac{{\lambda }^{2}+{a}^{2}+{b}^{2}}{ab}\;{x}_{1},\hfill \\ \hfill {x}_{2}& =\frac{{\lambda }^{2}+{a}^{2}}{ab}\;{x}_{1}.\hfill \end{align*}$$

The y₂ can be rewritten as ${y}_{2}={y}_{1}\left[{f}_{1}+{f}_{2}\right]$ which leads to a simpler form of all y_l's [40]

$$\begin{equation}{y}_{l}=\frac{\lambda }{a}\;{x}_{1}\;{T}_{l}.\end{equation} \tag{ C9 }$$

After a bit of algebra, one finds x_l to be

$$\begin{equation}{x}_{l}={x}_{1}\;\left[{T}_{l}-\frac{b}{a}{T}_{l-1}\right].\end{equation} \tag{ C10 }$$

So far we found the general solutions of the recursion formulas equations (C1)–(C4). The comparison of equations (C2) and (C3) leads to

$$\begin{equation}{x}_{\frac{N}{2}+1}=0,\end{equation} \tag{ C11 }$$

because the recursion formulas themselves do not care about any index limitation. The last equation means only that the wave function of a finite system has to vanish outside, at the boundary, yielding the quantization rule.

The extended states can be obtained with

$$\begin{align*}\hfill {f}_{1}& ={\mathrm{e}}^{2\text{i}\;kd},\hfill \\ \hfill {f}_{2}& ={\mathrm{e}}^{-2\text{i}\;kd},\hfill \end{align*}$$

where equations (C6) and (C7) relate kd and λ, and T_l is recast as

$$\begin{equation}{T}_{l}=\frac{\mathrm{sin}\left(2\;kd\;l\right)}{\mathrm{sin}\left(2\;kd\right)}.\end{equation} \tag{ C12 }$$

The last equation for T_l yields via equations (C10) and (C11) the quantization condition. Thus, the momenta k obey

$$\begin{equation*}\mathrm{tan}\left[kd\;\left(N+1\right)\right]=\frac{{\Delta}}{t}\;\mathrm{tan}\left(kd\right),\end{equation*}$$

where each solution defines two states with the energy E_± from equation (35).

The decaying states depend strongly on the interplay of Δ and t. The ansatz is

$$\begin{align*}\hfill {f}_{1}& =s\;{\mathrm{e}}^{2\;qd},\hfill \\ \hfill {f}_{2}& =s\;{\mathrm{e}}^{-2\;qd},\hfill \end{align*}$$

where s is defined as

$$\begin{equation*}s=\left\{\begin{matrix}\hfill +1,\hfill & \hfill \vert {\Delta}\vert { >}\vert t\vert \hfill \\ \hfill -1,\hfill & \hfill \vert t\vert { >}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation*}$$

Finally the coefficient T_l becomes

$$\begin{equation*}{T}_{l}\left(qd\right){:=}{s}^{l-1}\;\;\;\frac{\mathrm{sinh}\left(2\;qd\;l\right)}{\mathrm{sinh}\left(2\;qd\right)}.\end{equation*}$$

The proper q, if existent, leads to two states and satisfies

$$\begin{equation*}\mathrm{tanh}\left[qd\;\left(N+1\right)\right]=m\;\mathrm{tanh}\left(qd\right),\end{equation*}$$

where m is

$$\begin{equation}m\;{:=}\;\left\{\begin{matrix}\hfill \frac{{\Delta}}{t},\quad \;\mathrm{i}\mathrm{f}\;\vert {\Delta}\vert {\geqslant}\vert t\vert \hfill \\ \hfill \frac{t}{{\Delta}},\quad \;\mathrm{i}\mathrm{f}\;\vert t\vert {\geqslant}\vert {\Delta}\vert \hfill \end{matrix}\right..\end{equation} \tag{ C13 }$$

In total we have already all N non normalized states with respect to the chain α and this approach holds as long as |Δ| ≠ |t|.

The remaining cases start again from the equations (C1)–(C4).

Case 2. Eigenvectors at the Kitaev point. We consider now Δ = −t, or b = 0, and we have to solve

$$\begin{align*}\hfill a\;{y}_{1}& =\lambda \;{x}_{1},\hfill \\ \hfill -a\;{x}_{N/2}& =\lambda \;{y}_{N/2},\hfill \\ \hfill -a\;{x}_{l}& =\lambda \;{y}_{l},\hfill \\ \hfill a\;{y}_{l+1}& =\lambda \;{x}_{l+1},\hfill \end{align*}$$

where l runs from 1 to (N − 2)/2. A zero energy mode is obviously not existing on the α subchain, because λ = 0 would lead to ${\to {v}}_{\alpha }=0$ which is not an eigenvector by definition. These zero modes belong to the subchain β for Δ = −t. The only possible eigenvalues for the extended modes of the α chain are λ = ±2t [1, 2], see equation (35). Recalling a = −2it, leads to N/2 independent solutions of dimerised pairs (x_l, y_l) with y_l = ∓ix_l and the signs are with respect to the eigenvalues.

The last cases belong to Δ = t (a = 0), where we search for the solution of

$$\begin{align*}\hfill \lambda \;{x}_{1}& =0,\hfill \\ \hfill \lambda \;{y}_{N/2}& =0,\hfill \\ \hfill b\;{x}_{l+1}& =\lambda \;{y}_{l},\hfill \\ \hfill -b\;{y}_{l}& =\lambda \;{x}_{l+1},\hfill \end{align*}$$

where l runs from 1 to (N − 2)/2. The first (second) line clearly states that either λ is zero and/or x₁ (y_N/2). The zero λ means on the one hand that most entries vanish x₂ = x₃ = ⋯ = x_N/2 = 0 and y₁ = y₂ = ⋯ = y_(N−2)/2 = 0, since b = 2it ≠ 0 to avoid a trivial Hamiltonian. On the other hand we have two independent solutions, first

$$\begin{align*}\hfill {x}_{1}& =1,\hfill \\ \hfill {y}_{N/2}& =0,\hfill \end{align*}$$

and second

$$\begin{align*}\hfill {x}_{1}& =0,\hfill \\ \hfill {y}_{N/2}& =1,\hfill \end{align*}$$

describing the isolated MZM's at opposite ends of the chain. In the case of λ = ±2t, we have N − 2 independent solutions in form of pairs (y_l, x_l+1) with y_l = ±ix_l+1.

The non trivial solutions for ${\to {v}}_{\beta }$ follow by replacing ${x}_{l}\to {\mathcal{X}}_{l}$, ${y}_{l}\to {\mathcal{Y}}_{l}$ and t → −t everywhere.

The eigenvectors have similar shape

$$\begin{align*}\hfill {\to {v}}_{\alpha }& ={\left({x}_{1},\;{y}_{1},\;{x}_{2},\;{y}_{2},\dots ,\;{x}_{\frac{N-1}{2}},\;{y}_{\frac{N-1}{2}},\;{x}_{\frac{N+1}{2}}\right)}^{\text{T}},\hfill \\ \hfill {\to {v}}_{\beta }& ={\left({\mathcal{X}}_{1},\;{\mathcal{Y}}_{1},\;{\mathcal{X}}_{2},\;{\mathcal{Y}}_{2},\dots ,\;{\mathcal{X}}_{\frac{N-1}{2}},\;{\mathcal{Y}}_{\frac{N-1}{2}},\;{\mathcal{X}}_{\frac{N+1}{2}}\right)}^{\text{T}},\hfill \end{align*}$$

but the last entry is different compared to the even N case. Although both subchains have the same spectrum, it is possible to consider a superposition of eigenstates of the full Hamiltonian which belongs to only one chain, for example α. We consider ${\to {v}}_{\beta }$ to be zero.

The eigenvector system for ${\to {v}}_{\alpha }$ reads

$$\begin{align*}\hfill a\;{y}_{1}& =\lambda \;{x}_{1},\hfill \\ \hfill -b\;{y}_{\frac{N-1}{2}}& =\lambda \;{x}_{\frac{N+1}{2}},\hfill \end{align*}$$

and

$$\begin{align*}\hfill b\;{x}_{i+1}-a\;{x}_{i}& =\lambda \;{y}_{i},\hfill \\ \hfill a\;{y}_{l+1}-b\;{y}_{l}& =\lambda \;{x}_{l+1},\hfill \end{align*}$$

with $l=1,\;\;\dots ,\;\;\frac{N-3}{2}$ and $i=1,\;\;\dots ,\;\;\frac{N-1}{2}$.

If we consider a, b and λ all to be different from zero, we find again that the entries of ${\to {v}}_{\alpha }$ are Fibonacci polynomials obeying the same recursion formula as in the even N case and lead to the same solution

$$\begin{align*}\hfill {y}_{i}& =\frac{\lambda }{a}\;{T}_{i}\;{x}_{1},\hfill \\ \hfill {x}_{l}& =\left[{T}_{l}-\frac{b}{a}\;{T}_{l-1}\right]\;{x}_{1},\hfill \end{align*}$$

where $l=1,\;\;\dots ,\;\frac{N+1}{2}$ and T_{l, (i)} is as before. The ansatz f₁ = e^{2i kd}, f₂ = e^{−2i kd} for the extended states influences T_l (T_i analogously)

$$\begin{equation*}{T}_{l}=\frac{\mathrm{sin}\left(2\;kd\;l\right)}{\mathrm{sin}\left(2\;kd\right)},\end{equation*}$$

and leads via

$$\begin{equation*}{y}_{\frac{N+1}{2}}=0,\end{equation*}$$

to the equidistant quantization $k\;\equiv \;{k}_{n}=\frac{n\;\pi }{N+1}$ with n = 1, ..., (N − 1)/2, due to the number of eigenvectors of the single SSH-like chain. Both chains share the same spectrum for odd N and thus we have in total n = 1, ..., N, n ≠ (N + 1)/2.

We report here shortly on all other parameter situations.

• (a)  
  If we consider a and b to be different from zero, but λ = 0, we find only one state

  $$\begin{equation}{x}_{l+1}={\left(\frac{{\Delta}-t}{{\Delta}+t}\right)}^{l}\;{x}_{1},\end{equation} \tag{ C14 }$$

  and l runs from 1 to (N − 1)/2.
• (b)  
  If Δ = t, i.e. a = 0, but λ = ±2t ≠ 0, we find (N − 1)/2 solutions (y_l, x_l+1) with y_l = ±ix_l+1, l = 1, ..., (N − 1)/2 and x₁ = 0 for all.The zero mode of this setting (Δ = t) is an MZM localized on x₁ = 1, while all other components are zero.
• (c)  
  If Δ = −t (b = 0) and λ ≠ 0 we find (N − 1)/2 solutions of the form (x_l, y_l) with y_l = ±ix_l l = 1, ..., (N − 1)/2 and ${x}_{\frac{N+1}{2}}=0$ for all of them. The MZM is localised at ${x}_{\frac{N+1}{2}}=1$ for b = 0.The results for the α chain follow again by replacing ${x}_{l}\to {\mathcal{X}}_{l}$, ${y}_{l}\to {\mathcal{Y}}_{l}$ and t → −t.

The most general form for μ is

$$\begin{equation}\mu =2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;f\left(\theta \right),\end{equation} \tag{ F9 }$$

where the function f(θ) accounts for all possible ratios of μ and $\sqrt{{t}^{2}-{{\Delta}}^{2}}$. The case (a) enforces the function f(θ) to be real valued, because both μ and $\sqrt{{t}^{2}-{{\Delta}}^{2}}$ are real. Further, we find that

$$\begin{equation}{f}^{2}\left(\theta \right){\leqslant}1,\end{equation} \tag{ F10 }$$

since 4(t² − Δ²) − μ² ⩾ 0. Please note that equation (F10) needs only to hold for θ on a finite set. From all possible functions f(θ), a convenient choice is f = cos(θ). The reason behind our specific choice is the form of R_1,2, because f leads in $\sqrt{4\;ab\;-{\mu }^{2}}$ to

$$\begin{align*}\hfill \sqrt{4\;ab\;-{\mu }^{2}}\;& =\;\sqrt{4\;\left({t}^{2}-{{\Delta}}^{2}\right)\;-{\mu }^{2}}\hfill \\ \hfill & =\;\sqrt{4\;\left({t}^{2}-{{\Delta}}^{2}\right)\;\left[1-{\mathrm{cos}}^{2}\left(\theta \right)\right]}\hfill \\ \hfill & =\;2\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{sin}\left(\theta \right),\hfill \end{align*}$$

and R_1,2(f) become

$$\begin{align*}\hfill {R}_{1,2}\left(f\right)\;& =\;\frac{-i\mu {\pm}\sqrt{4\;ab\;-{\mu }^{2}}}{2}\hfill \\ \hfill & =\;\frac{-2i\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{cos}\left(\theta \right){\pm}2\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{sin}\left(\theta \right)}{2}.\hfill \end{align*}$$

Simplifications lead to

$$\begin{equation*}{R}_{1,2}\left(f\right)=-i\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\left\{\begin{matrix}\hfill {\mathrm{e}}^{\text{i}\theta },\hfill \\ \hfill {\mathrm{e}}^{-\text{i}\theta }\hfill \end{matrix}\right..\end{equation*}$$

Let us focus on the determinant. We find ${R}_{1}^{j}-{R}_{2}^{j}$ to be

$$\begin{equation*}{R}_{1}^{j}-{R}_{2}^{j}={\left[-i\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\right]}^{j}\;2i\;\mathrm{sin}\left(\theta \;j\right).\end{equation*}$$

Consequently the determinant reads

$$\begin{equation}\mathrm{det}\left({\mathcal{H}}_{c}\right)={\left({t}^{2}-{{\Delta}}^{2}\right)}^{N}\;{\left[\frac{\mathrm{sin}\left[\theta \;\left(N+1\right)\right]}{\mathrm{sin}\left(\theta \;\right)}\right]}^{2},\end{equation} \tag{ F11 }$$

and vanishes for θ = nπ/(N + 1) (n = 1, ..., N) or t² = Δ². Since Δ, t and θ define together with f₁ = cos(θ) the chemical potential, we find that the determinant of the Kitaev chain is zero if, and only if:

• (a)  
  $\mu =2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{cos}\left(\frac{n\;\pi }{N+1}\right)$,
• (b)  
  μ = 0 and t² = Δ²,

for n = 1, ..., N, t² ⩾ Δ² and for all N. A feature of odd N is the value n = N + 1/2 yielding θ = π/2, i.e. μ = 0 for all values of Δ, t for t² ⩾ Δ². In fact μ = 0 holds for odd N everywhere, as we already know from previous discussion in appendix B.

We found all zeros in case (a) and we continue with (b).

We follow the same way of argumentation as above, but we have to keep in mind that t² − Δ² ⩾ 0, and 4(t² − Δ²) < μ². The first step is to reshape the square root in R_1,2

$$\begin{equation}\sqrt{4\;ab-{\mu }^{2}}=i\sqrt{{\mu }^{2}-4\;\left({t}^{2}-{{\Delta}}^{2}\right)},\end{equation} \tag{ F12 }$$

where we find a similar situation as in the previous scenario. Our ansatz is

$$\begin{equation}\mu =2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}g\left(\theta \right),\end{equation} \tag{ F13 }$$

where the function g(θ) is real and obeys

$$\begin{equation}{g}^{2}\left(\theta \right){\geqslant}1,\end{equation} \tag{ F14 }$$

since μ² ⩾ 4(t² − Δ²). The candidates of our choice are g_±(θ) = ±cosh(θ), where θ is real. The square root becomes now

$$\begin{equation*}\sqrt{{\mu }^{2}-4\;\left({t}^{2}-{{\Delta}}^{2}\right)}=2\;\sqrt{\left({t}^{2}-{{\Delta}}^{2}\right)}\;\mathrm{sinh}\left(\theta \right),\end{equation*}$$

and we find R_1,2(g₊) to be

$$\begin{align*}\hfill {R}_{1,2}\left({g}_{+}\right)\;& =\;\frac{-i\mu {\pm}i\sqrt{{\mu }^{2}-4\;\left({t}^{2}-{{\Delta}}^{2}\right)}}{2}\hfill \\ \hfill & =\;-i\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\left[\mathrm{cosh}\left(\theta \right)\mp \mathrm{sinh}\left(\theta \right)\right].\hfill \end{align*}$$

Simplifications yield

$$\begin{equation*}{R}_{1,2}\left({g}_{+}\right)=-i\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\left\{\begin{matrix}\hfill {\mathrm{e}}^{-\theta }\hfill \\ \hfill {\mathrm{e}}^{\;\theta }\hfill \end{matrix}\right.,\end{equation*}$$

and the determinant becomes:

$$\begin{equation}\mathrm{det}\left({\mathcal{H}}_{c}\right)={\left({t}^{2}-{{\Delta}}^{2}\right)}^{N}\;{\left[\frac{\mathrm{sinh}\left(\theta \left[N+1\right)\right]}{\mathrm{sinh}\left(\theta \right)}\right]}^{2}.\end{equation} \tag{ F15 }$$

The determinant vanishes only if t² = Δ², which by virtue of equation (F13) implies μ = 0, because the fraction of the hyperbolic sine functions is always positive. The use of g_1,− = −cosh(θ) leads to equation (F15) again.

We consider here Δ² ⩾ t² and 4(t² − Δ²) − μ² ⩽ 0. We start by manipulating the square root in R_1,2

$$\begin{equation}\sqrt{4\;ab-{\mu }^{2}}=\sqrt{4\;\left({t}^{2}-{{\Delta}}^{2}\right)-{\mu }^{2}}=i\;\sqrt{{\mu }^{2}+4\;\left({{\Delta}}^{2}-{t}^{2}\right)}.\end{equation} \tag{ F16 }$$

Our ansatz is $\mu =2\;\sqrt{{{\Delta}}^{2}-{t}^{2}}\;v\left(\theta \right)$ with a real valued function v(θ), without further restrictions, because

$$\begin{equation*}{\mu }^{2}=4\;\left({{\Delta}}^{2}-{t}^{2}\right)\;{v}^{2}\left(\theta \right){\geqslant}-4\;\left({{\Delta}}^{2}-{t}^{2}\right),\end{equation*}$$

in view of μ² ⩾ −4(Δ² − t²). The square root in R_1,2 becomes in general

$$\begin{equation*}i\;\sqrt{{\mu }^{2}+4\;\left({{\Delta}}^{2}-{t}^{2}\right)}=i\;\left({{\Delta}}^{2}-{t}^{2}\right)\;\sqrt{{v}^{2}\left(\theta \right)+1},\end{equation*}$$

and one sees immediately that v(θ) = sinh(θ), $\theta \in \mathbb{R}$ is an appropriate choice. We find for R_1,2 the form

$$\begin{equation*}{R}_{1,2}\left(v\right)=-i\;\sqrt{{{\Delta}}^{2}-{t}^{2}}\;\left\{\begin{matrix}\hfill -{\mathrm{e}}^{-\theta }\hfill \\ \hfill {\mathrm{e}}^{\theta }\hfill \end{matrix}\right.,\end{equation*}$$

where the negative sign in front of the exponential forces us to distinguish between even and odd N. The determinant reads finally

$$\begin{equation*}\mathrm{det}\left({\mathcal{H}}_{c}\right)={\left({{\Delta}}^{2}-{t}^{2}\right)}^{N}\;\left\{\begin{matrix}\hfill \frac{\mathrm{cosh}\left[\theta \left(N+1\right)\right]}{\mathrm{cosh}\left(\theta \right)},\hfill & \hfill N\;\text{even}\hfill \\ \hfill \frac{\mathrm{sinh}\left[\theta \left(N+1\right)\right]}{\mathrm{sinh}\left(\theta \right)},\hfill & \hfill N\;\text{odd}\hfill \end{matrix}\right.,\end{equation*}$$

and it is never zero, except for Δ² = t² at μ = 0.

In summary, for μ ≠ 0, we have only non trivial, zero determinants in case (a). How can one be sure that no zero is missed especially in the settings (b) and (c)? This follows immediately from equation (F8), because the determinant vanishes only if

$$\begin{equation*}{R}_{1}^{N+1}={R}_{2}^{N+1}.\end{equation*}$$

Consequently we need first of all |R₁| = |R₂|. The second part is to find the proper phase factors and all of them lie on a circle with radius |R₁| in the complex plane. We have found non trivial solutions only for scenario (a).

In total, we found all conditions $\mathrm{det}\left({\mathcal{H}}_{\text{KC}}\right)=0$. The general case is when the chemical potential is

$$\begin{equation}\mu =2\;\sqrt{{t}^{2}-{{\Delta}}^{2}}\;\mathrm{cos}\left(\frac{n\;\pi }{N+1}\right),\end{equation} \tag{ F17 }$$

with t² ⩾ Δ² and n = 1, ..., N, i.e. the chemical potential corresponds to the energy levels of a linear chain with hopping $\sqrt{{t}^{2}-{{\Delta}}^{2}}$. The case μ = 0 and t² = Δ² is included in equation (F17).

Further, the determinant of a Kitaev chain with odd number of sites is zero if μ = 0 for all values of Δ and t.