Extended Kitaev chain with longer-range hopping and pairing

TL;DR

This paper explores an extended Kitaev chain model with long-range interactions, revealing complex topological phases and multiple Majorana modes at edges, using topological invariants, exact diagonalization, and transfer matrix methods.

Contribution

It introduces a generalized transfer matrix approach and analyzes the effects of long-range hopping and pairing on Majorana modes and topological phases.

Findings

Multiple Majorana modes can appear at each edge with extended couplings.

Rich phase diagrams depend on the range of interactions and symmetry considerations.

Generalized transfer matrix method effectively calculates edge Majorana modes.

Abstract

We consider the Kitaev chain model with finite and infinite range in the hopping and pairing parameters, looking in particular at the appearance of Majorana zero energy modes and massive edge modes. We study the system both in the presence and in the absence of time reversal symmetry, by means of topological invariants and exact diagonalization, disclosing very rich phase diagrams. In particular, for extended hopping and pairing terms, we can get as many Majorana modes at each end of the chain as the neighbors involved in the couplings. Finally we generalize the transfer matrix approach useful to calculate the zero-energy Majorana modes at the edges for a generic number of coupled neighbors.

Equations

$$H$$

$$w_{\ell }=w_0{d}_{\ell }^{-\alpha },{\Delta }_{\ell }=\Delta d_{\ell }^{-\beta },$$

$$H$$

$${\mathcal{C}}^{-1}\mathcal{H}(k)\mathcal{C}=-\mathcal{H}(-k).$$

$$\text{h}(k)=\left(0,h_y(k),h_z(k)\right)$$

$$\mathcal{H}(k)=\text{h}(k)\cdot \mathbf{\sigma }$$

$$\mathsf{w}=\frac{1}{2\pi }\oint d{\theta }_k=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }dk\frac{{\partial }_k{\hat{h}}_z(k)}{{\hat{h}}_y(k)}$$

$$\upsilon =\text{sign}(h_z(0)h_z(\pi ))$$

$$-2w_0{H}_r^{\alpha }<\mu <2w_0{R}_r^{\alpha }$$

$$H_r^{\alpha }=\sum _{\ell =1}^{r}1/{\ell }^{\alpha }$$

$$R_r^{\alpha }=\sum _{\ell =1}^r(-1{)}^{\ell +1}/{\ell }^{\alpha }$$

$$h_y(k)$$

$$h_z(k)$$

$$k_{2n}=\frac{2\pi n}{r},k_{2n+1}=\frac{2\pi (n+1)}{r+1}.$$

$$h_z(0)=-\mu /2-r{w}_0$$

$$h_z(k_{2n})=-\mu /2$$

$$h_z(k_{2n+1})=-\mu /2-w_0$$

$$\mathsf{w}=1,$$

$$\mathsf{w}=r,$$

$$-2\sum _{\ell =1}^r{w}_{\ell }\cos {\phi }_{\ell }<\mu <-2\sum _{\ell =1}^r(-1{)}^{\ell }{w}_{\ell }\cos {\phi }_{\ell }$$

$$-2\sum _{\ell =1}^r\frac{\cos {\phi }_{\ell }}{{\ell }^{\alpha }}<\frac{\mu }{w_0}<-2\sum _{\ell =1}^r(-1{)}^{\ell }\frac{\cos {\phi }_{\ell }}{{\ell }^{\alpha }}.$$

$$E_k^{\pm }=\sum _{\ell =1}^r{w}_{\ell }\sin {\phi }_{\ell }\sin (k\ell )$$

$$\displaystyle\pm\sqrt{\Big{[}\sum_{\ell=1}^{r}w_{\ell}\cos\varphi_{\ell}\cos(k\ell)+\frac{\mu}{2}\Big{]}^{2}+\Big{[}\sum_{\ell=1}^{r}\Delta_{\ell}\sin(k\ell)\Big{]}^{2}}$$

$$\eta \equiv \underset{\{k\}}{\min }(E_k^{+}{E}_k^{-})\underset{\{k\}}{\max }(E_k^{+}{E}_k^{-})$$

$$\mathcal{H}(k)=\left(\begin{array}{cc}{h}_0(k)+h_z(k)& -i{h}_y(k)\\ i{h}_y(k)& h_0(k)-h_z(k)\end{array}\right)$$

$$h_0(k)\pm h_z(k)=$$

$$\mp \frac{1}{2}\left\{\mu +w_0\left[e^{-i{\phi }_0}{\text{Li}}_{\alpha }(e^{-i(\phi \pm k)})+e^{i{\phi }_0}{\text{Li}}_{\alpha }(e^{i(\phi \pm k)})\right]\right\}$$

$$h_y(k)=i\frac{\Delta }{2}\left[{\text{Li}}_{\alpha }(e^{-ik})-{\text{Li}}_{\alpha }(e^{ik})\right]$$

$${\mu }_{c1}=-2w_0\zeta (\alpha )$$

$${\mu }_{c2}=2w_0\eta (\alpha )=2w_0(1-2^{1-\alpha })\zeta (\alpha )$$

$$h_z(k)=h_z(2\pi -k),$$

$$h_y(k)=-h_y(2\pi -k).$$

$${\mu }_{c1}=-w_0\left[{\text{Li}}_{\alpha }(e^{-i{\phi }_0})+{\text{Li}}_{\alpha }(e^{i{\phi }_0})\right]$$

$${\mu }_{c2}=-w_0\left[{\text{Li}}_{\alpha }(-e^{-i{\phi }_0})+{\text{Li}}_{\alpha }(-e^{i{\phi }_0})\right]$$

$$c_{2j-1}=a_{j}+a^{\dagger}_{j},\;\;\;\;\;c_{2j}=i\big{(}a^{\dagger}_{j}-a_{j}\big{)}$$

$$H$$

Full text

I Introduction

Topological orders and their related quantum phenomena have attracted a lot of attention in the last few years Wen . One of the most interesting model which shows a topological order is the Kitaev chain Kitaev . This model describes a p-wave superconducting wire that, under certain conditions, has a gapped bulk together with zero energy unpaired Majorana modes localized at the edges of the system, which are robust against disorder Chen , local impurities Ising or dynamical perturbations DeGottardi2 ; Zvyagin . The presence of such modes has recently been observed experimentally ExperimentalMFs ; Exp2 and optical implementations of generalized Kitaev model have been proposed optics . This one-dimensional system represents the simpliest playground for designing topological quantum computers Alicea ; Akhmerov , because of its straightforward mapping to a spin system Ising ; Pfeuty ; Lieb .

An extension of the Kiteav model has been recently proposed by considering a longer range in the hopping and pairing terms appearing in the Hamiltonian Ising ; DeGottardi ; Ghazaryan ; Lepori ; Lepori2 ; delgado ; Lepori3 ; Shi ; Pientkaprb ; Pientka ; pachos . We will consider such a situation for both long-range and finite-range interactions. In the latter case, if time reversal symmetry (TRS) holds, one can get many Majorana zero modes (MZMs) per edge as in the case of multiband spin-orbit coupled superconductor nanowire with parallel Zeeman splitting Tewari . In this case, according to standard symmetry classification classification , the system belongs to the class BDI, and the number of zero modes is dictated by the $\mathbb{Z}$-valued topological invariant (TI) winding number ${\sf w}$ Ising ; Zinvariant1 . Interestingly, for certain range of parameters, the number of Majorana modes can be equal to the number of neighbors. If time reversal symmetry is broken, the system belongs to class $D$ and one can have at most one Majorana zero mode per edge Ising in the regime of parameters where the pfaffian $\mathbb{Z}_{2}$ topological invariant Kitaev is not trivial. It has been shown that such invariant, related to the fermionic parity of the ground state, corresponds to the quantized Zak-Berry phase Budich .

In this paper we study in a systematic way the phase diagrams of the extended Kitaev model, considering also the long-range limit both in the presence and in the absence of time reversal symmetry, by calculating the topological invariants and performing exact diagonalization of the Hamiltonian. We also provide a detailed derivation of the transfer matrix approach in the case where a generic number of neighbors are involved in the couplings and many Majorana states are located at the edges, consistently with the topological invariant analysis.

II Long-Range Hamiltonian

We propose an extended Kitaev chain model taking into account $r$ neighbor interactions in the hopping and pairing terms. We assume both of these interactions to be algebraically decreasing with the distance between two different lattice sites. The model with infinitely long-range pairing has been studied in Ref. Lepori . This model shows, for certain physical regimes, Majorana zero modes and massive edge modes identified as topological massive Dirac fermions in Ref. delgado , because their energies are separated from the excited states by a finite gap even in the thermodynamic limit.

The fermionic Hamiltonian we will consider, which generalizes the Kitaev chain, is the following

[TABLE]

where $L$ is the number of lattice sites and $\mu$ the chemical potential. The extended hopping and pairing coupling terms, $w_{\ell}$ and $\Delta_{\ell}$, can be generic variables, although, in what follows, we will assume the following form for those parameters

[TABLE]

which couple the lattice site $j$ with the site $(j+\ell)$. The hopping parameter $w_{\ell}$ can acquire a phase, $w_{\ell}e^{i\varphi_{\ell}}$, with $\varphi_{\ell}$, $w_{0}$, $\Delta$ real values, in the case of broken time reversal symmetry. The index $\ell$ runs over the neighbor sites and $d_{\ell}$ is $d_{\ell}=\min(\ell,L-\ell)$ for closed boundary conditions and $d_{\ell}=\ell$ for open boundary conditions. The exponents $\alpha$ and $\beta$ characterize the rate of decay for the parameters if they are assumed not negative. For $\alpha\rightarrow\infty$, $\beta\rightarrow\infty$ one recovers the standard short-range Kitaev model. Let us suppose now to close the chain with periodic (PBC) or antiperiodic boundary conditions (ABC) and make the Fourier transform $a_{j}=\frac{1}{\sqrt{L}}\sum_{k}a_{k}e^{-ikj}$ with $k=\frac{2\pi n}{L}$ for PBC and $k=\frac{2\pi n+\pi}{L}$ for ABC, with $n$ integer numbers. The Hamiltonian can be written as follows

[TABLE]

The spectrum has a particle-hole symmetry since, under particle-hole transformation, ${\cal C}=\sigma_{x}{\cal K}$ (where ${\cal K}$ is the complex conjugation operator),

[TABLE]

The time reversal condition ${\cal H}(-k)^{*}={\cal H}(k)$ is satisfied only if we consider real hopping terms $w_{\ell}$, thus if $\varphi_{\ell}=0,n\pi$. In the next sections we will analyze both situations, with and without TRS.

Before we proceed a comment about boundary conditions is in order. In Ref. Lepori , ABC are assumed to preserve pairing terms in $H$, on the same time this choice would destroy long-range hopping terms for finite $L$ and $r\geq L/2$. For PBC the reverse holds. If $r<L/2$, instead, both choices preserve long-range interaction terms and translational invariance. Closed chains with ABC and PBC, however, may give different results. For finite $r$, in the expression for ${\cal H}(k)$, terms like $\cos\left(\ell\left(\frac{2\pi n}{L}+\frac{\pi}{L}\right)\right)$ and $\sin\left(\ell\left(\frac{2\pi n}{L}+\frac{\pi}{L}\right)\right)$ for ABC can be confused with the corresponding PBC conterparts, $\cos\left(\ell\left(\frac{2\pi n}{L}\right)\right)$ and $\sin\left(\ell\left(\frac{2\pi n}{L}\right)\right)$, in the limit $L\rightarrow\infty$, if $\ell$ remains finite. If, instead, we consider an infinite number of interacting neighbors (always in the thermodynamic limit, $L\rightarrow\infty$), $d_{\ell}$ is always $\ell$ and the above terms in (3) generates polylogarithmic functions $Li_{\alpha}(e^{\pm ik\ell})$ where $k$ becames a continuum variable. In the thermodynamic limit, therefore, ${\cal H}(k)$ is always the same, for both PBC and ABC.

III Topological invariants

Topological phases are described by different topological invariants according to the presence or the absence of time reversal symmetry. If time reversal is preserved ($\varphi_{\ell}=0$) there is a $\mathbb{Z}$ topological invariant given by the winding number. In this case we can define the vector

[TABLE]

such that ${\cal H}(k)$ in Eq. (3), for $\varphi_{\ell}=0$, can be written as

[TABLE]

where $\bm{\sigma}$ is the vector made of Pauli matrices. In terms of unit vector $\hat{\textbf{h}}(k)=\textbf{h}(k)/|\textbf{h}(k)|$ we define the winding number

[TABLE]

In the case of broken TRS ($\varphi_{\ell}\neq 0$), as we can see from Eq. (3), we also have to consider the component of ${\cal H}(k)$ proportional to $\mathds{1}$ (or by a transformation, proportional to ${\sigma}_{x}$) so that the winding number ${\sf w}$ cannot be defined. For $k=0,\pm\pi$ we have only the Hamiltonian term along ${\sigma}_{z}$, while the term along ${\sigma}_{y}$ vanishes. In such a situation a good topological invariant is given by

[TABLE]

which can be used also for time reversal symmetry cases.

IV Finite number of neighbors: topological phase diagrams

Here we will focus on the case of finite number of interacting neighbors. At first we will allow the hopping term to be longer ranged, keeping the pairing interaction short-ranged, and then we will consider the opposite case. At the end we will allow both terms to be longer ranged, analyzing the topological phase in the presence or in the absence of TRS.

IV.1 Extended-range hopping

We will now analyze Eq. (II) in the limit of only extended-range hopping, i.e. $\beta\rightarrow\infty$, for finite $r$. In Fig. 1 we observe that the topological regime (${\sf w}=\pm 1$) of parameters increases with $\alpha$, the exponent of the hopping term. The result we found for this specific Hamiltonian is consistent with the one obtained in delgado . Only one MZM per edge can be found, since the dependence of the pairing term on $k$ goes as $\sin(k)$ and not as $\sin(\ell k)$ ($\ell=1,\dots,r$). This means that moving along the first Brillouin zone (in $k$-space), the winding vector makes only one cycle around $(0,0)$. To better show this aspect we report in Fig. 2(b) the points spanned by the unormalized winding vector $\textbf{h}(k)=\left(0,\Delta\sin(k),-\mu/2-w_{0}\sum_{\ell}\ell^{-\beta}\cos(\ell k)\right)$. At any finite number of neighbors $r\leq L/2$ and exponent $\alpha$, the topological regime for $w_{0}>0$ is defined by the condition

[TABLE]

where

[TABLE]

Here $H_{r}^{\alpha}$ is the generalized harmonic number of order $\alpha$ of $r$, and $R_{r}^{\alpha}=(1-2^{1-\alpha})\zeta(\alpha)-\zeta(\alpha,1+r)+2^{1-\alpha}\zeta(\alpha,1+\lfloor r/2\rfloor)$, with $\zeta(\alpha,n)=\sum_{\ell=0}^{\infty}1/(\ell+n)^{\alpha}$ the generalized Riemann zeta function. For $\alpha\rightarrow\infty$ (short-range limit), $H_{r}^{\infty}=R_{r}^{\infty}=1$, and Eq. (9) reduces to the standard Kitaev condition. In the case of infinite number of neighbors ($r\rightarrow\infty$), in the thermodynamic limit, we get $H_{\infty}^{\alpha}=\zeta(\alpha)$, the Riemann zeta function, and $R_{\infty}^{\alpha}=\eta(\alpha)=(1-2^{1-a})\zeta(\alpha)$ the Dirichlet eta function.

IV.2 Extended-range pairing

We now analyze longer range Hamiltonian with a finite number $r$ of neighbors only in the pairing term, namely $\alpha\rightarrow\infty$ and choosing a finite $\beta$. Again the winding number, as before, can take only the values $0,\pm 1$, however, we observe an alternation, inside the topological phase, of the values $\pm 1$ (see Fig. 3). This behavior is related to the nodes along h_{y}(k)=\Delta\sum_{\ell=1}^{r}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sin(k\ell)/{\ell}^{\beta}}=0, as shown in Fig. 2(c), so that the origin can be surrounded clock- or counterclockwise. Notice that the full size of the topological phase is the same as that of the short-range model. This behavior, not present in the extended hopping Hamiltonian, is more emphatic as $r$ increases (we are always assuming $r$ to be finite) while it disappears as $\beta\rightarrow\infty$, in which case we recover the topological phase diagram of the usual Kitaev model Kitaev .

IV.3 Extended-range hopping and pairing, with and without TR symmetry

We will look now to the topological phases when considering longer range in the hopping and in the pairing at the same time, assuming $\beta=\alpha$, in the presence or in the absence of time reversal symmetry.

IV.3.1 With time reversal symmetry

After considering separately the extended range in the hopping and in the pairing, we now include them together. The winding numbers which characterize the topological phases provide a way to count the edge modes of the chain in the open configuration readegreen . In the special case where both the hopping and the pairing terms involve many neighbors ($r\geq 1$) we observe that the winding number, and therefore, the number of Majorana modes at each edge, can be larger than one, up to $r$. In order to prove this result, let us consider uniform case ($\alpha=\beta=0$) of $r$-neighbor hopping and pairing. In this case we have

[TABLE]

The $2r$ zeros of $h_{y}(k)$, namely $k_{n}$ such that $h_{y}(k_{n})=0$, can be ordered, $k_{0}<k_{1}<k_{2}<k_{3}<....<k_{2r-1}$, and are

[TABLE]

The corresponding values of $h_{z}(k_{n})$ are

[TABLE]

In the regime of parameters where $h_{z}(k_{n})$ has alternate signs for ordered $\{k_{n}\}$, varying $k$ from [math] to $2\pi$ the winding vector surrounds the origin $r$ times. As a result, for $\alpha=\beta=0$, $w_{0}>0$, we get from Eq. (IV.3.1) the following topological phases

[TABLE]

and ${\sf w}=0$ otherwise. This result is verified numerically in Fig. 4 a) for $r=2$ where we obtain $\left|{\sf w}_{max}\right|=2$ which should correspond to having $2$ MZM per edge readegreen . For a generic value of $\alpha$ the calculation is more involved, but we observe that by increasing $\alpha$, the range of parameters in which ${\sf w}$ is maximum decreases and for $\alpha\rightarrow\infty$ we recover the phase diagram of the standard Kitaev chain model. However, increasing $\alpha$, in some range of parameters, ${\sf w}$ can take intermediate values between $1$ and $r$ with steps of $2$, namely if $r$ is even ${\sf w}$ can be equal to [math] (trivial phase), $1$ or any even number $\leq r$, while if $r$ is odd, ${\sf w}$ can take the values [math], $1$ or any odd number $\leq r$ (see Fig. 5, where we parametrize $\mu/\Delta=\cos\gamma$ and $w_{0}/\Delta=\sin\gamma$ and plot ${\sf w}$ varying $\gamma$ for $\alpha=\beta=0.05$ and different values of $r$). This behavior is due to the mirror symmetry of the map $h_{z}(h_{y})$ with respect to $h_{y}=0$, as one can see in Fig. 2(d). This means that the edge modes are created or annihilated in pair at each edge, with the only exception of the last Majorana modes. This behavior can be explained by the fact that increasing $\alpha$ for some range of parameters the high degeneracy at zero energy level is hardly sustained and it is easier for two Majorana modes at the same edge to annihilate each other because of the overlap. Since the system is particle-hole symmetric the same process occurs also on the other edge. If $r$ is even the last four zero modes can annihilate each other on each edge (going to the trivial phase) or by one per edge, by anti-dimerization, inceasing $\mu$. The last two Majorana modes, one at each edge, when ${\sf w}=1$, instead, are more robust because their annihilation requires an overlap between the two wavefunctions peaked at long (infinite) distance as in the usual Kitaev chain.

IV.3.2 With broken time reversal symmetry

Including broken time reversal effects implies coupling between odd as well as even index Majorana operators so that only topological phase with an odd ${\sf w}$ will survive because of particle-hole symmetry of the spectrum DeGottardi . The time reversal symmetry is broken by the phase $\varphi_{\ell}$ in Eq. (3). We are going to study the effects of two different forms of this parameter: i) $\varphi_{\ell}=\varphi_{0}$, $\forall\ell=1,\dots,r$ and ii) $\varphi_{\ell}=\varphi_{0}\ell$ with $\ell=1,\dots,r$. Here we use the $\mathbb{Z}_{2}$ topological invariant $\upsilon=\operatorname{sign}\left(h_{z}(0)h_{z}(\pi)\right)$ to detect the topological phase, where $h_{z}$ is the term multplying $\sigma_{z}$ in Eq. (3). It has two possible values, $\upsilon=-1$ if the phase is topological and $\upsilon=1$ if it is trivial. The topological phase is permitted quite in general only if the condition

[TABLE]

is fulfilled, which, in our case, becomes

[TABLE]

For $\alpha=\beta=0$ and $\varphi_{\ell}=\varphi_{0}$ (if the phase does not depend on $\ell$, as in the plots on the r.h.s. of Fig. 6), then Eq. (19) is given by 2r\cos\varphi_{0}<\mu/w_{0}<2\cos\varphi_{0}\big{(}1+(-1)^{r}\big{)}/2, while for $\alpha=\beta=0$ and $\varphi_{\ell}=\varphi_{0}\,\ell$ (as in all the plots shown for different $\varphi_{0}$ on the l.h.s. of Fig. 6) then Eq. (19) becomes -2\cos\big{[}(r+1)\varphi_{0}/2\big{]}\sin(r\varphi_{0}/2)/\sin(\varphi_{0}/2)<\mu/w_{0}<\left\{1+\cos\big{[}(r+1)(\varphi_{0}+\pi)\big{]}+\tan(\frac{\varphi_{0}}{2})\sin\big{[}(r+1)(\varphi_{0}+\pi)\big{]}\right\}. The parametric regime where topological states is permitted shrinks as $\varphi_{0}$ increases, and disappears completely for $\varphi_{0}=\pi/2$, as one can easily check by putting this value on the expressions above.

Let us now write the bulk energy-momentum dispersion

[TABLE]

We find extended two-dimensional critical regions of parameters where $E_{k}^{\pm}$ vanishes for some values of $k$ (white regions in Fig. 6), as in the case of short-range Kitaev chain with broken time reversal symmetry DeGottardi . In particular, for $\varphi_{\ell}=\varphi_{0}\ell$, we can have also disconnected critical regions. These quantum critical points can be defined as the values of $\mu/w_{0}$ and $\Delta/w_{0}$ such that the following quantity

[TABLE]

is negative. In other words the critical regions can be defined by $\operatorname{sign}(\eta)=-1$, while the system is gapped if $\operatorname{sign}(\eta)=1$.

V Infinite number of neighbors: topological phase diagrams

Here we will consider an infinite number of neighbors, $r\rightarrow\infty$, first only in the pairing term, as done in Refs. Lepori2 and delgado , extending the analysis also for broken time reversal symmetry, and, afterwards, considering long range in both hopping and pairing terms. We will consider open boundary conditions, and perform exact diagonalization of the Hamiltonian Eq. (II), with $r,\,L\rightarrow\infty$ (numerically $L\gg 1$). Looking at the lowest energy levels, we find that, together with massless Majorana modes, we obtain massive edge states, as predicted for long-range pairing Lepori2 , separated by a finite gap, therefore also called topological Dirac fermions delgado .

V.1 Long-range pairing, with and without TR symmetry

In this section we will reconsider the Kitaev chain with long-range pairing, $w_{\ell}=w_{0}\delta_{\ell,1}$ and $\Delta_{\ell}=\Delta d_{\ell}^{-\alpha}$, in the presence of time reversal symmetry ($\varphi_{\ell}=\varphi_{0}=0$), or in its absence ($\varphi_{0}\neq 0$)

V.1.1 With time reversal symmetry

Let us distinguish three regimes of values for $\alpha$, which correspond to the regimes where i) ${\cal H}$ and its derivates are not defined in $k=0$, ii) only the derivates of ${\cal H}$ are not defined in $k=0$, and iii) both the Hamiltonian and its derivates are defined over all the Brilloiun zone. For $\alpha<1$ due to the discontinuity in $k=0$ of ${\cal H}(k)$, the winding number ${\sf w}$ as defined in Eq. (7), takes semi-integer values ${\sf w}=\pm\frac{1}{2}$ (see more comments about such winding numbers in Sec. V.2.1). For $\mu<2$ and $\alpha<1$ (here we choose $w_{0}=1$), the system has massive Dirac fermions at the edges, called massive edge modes (MEM), topologically protected by fermion parity and by a finite gap from bulk excitations, as shown in Ref. delgado . The topological nature of the non-local massive Dirac fermions purely comes from the long-range deformation of the original Kitaev model. For $1<\alpha<1.5$ there is a cohexistence of massless edge modes (Majorana zero modes, MZM) and MEM, as already pointed out in Ref. delgado . In this crossover region, we would like to give a more precise collocation of such massive and massless edge modes in the parameter space. We note that the presence of MEM also in the regions with ${\sf w}=0$ and ${\sf w}=1$, for points closed to $(\mu,\alpha)=(-2,1)$, as shown by the red-triangular points in Fig. 7, are consistent with the numerical scaling analysis for the mass (ground state energy $\Lambda_{0}\equiv E_{1}$) and the first gap ($\Delta E_{21}=E_{2}-E_{1}$) (see Fig. 8 for some examples). Quite interestingly, we find that in many cases there are finite gaps in the thermodynamic limit also between a few of energy levels, see Figs. 8 (b), (d). We also checked other points (dark-yellow square points in Fig. 8) where the mass is quite small, compatible with MZM. In the region signed by a blue star in Fig. 8 the scaling analysis is more difficult and not conclusive, therefore we cannot say anything about the edge modes there. However, this finding suggests that the winding number seems to be not enough to detect, for infinite-range model, the appearance of MZM, since at ${\sf w}=1$ one can find also MEM. For $\mu>2$ and $\alpha<1$ (blue region in Fig. 8) we get ${\sf w}=-1/2$ which seems therefore topologically not equivalent to the trivial phase with ${\sf w}=0$, for $\alpha>1$, (see comments on this in Sec. V.2.1) even if in both cases there are no edge modes (no EM). For $\alpha>3/2$ the system behaves like the short-range Kitaev model, with MZM for $-2<\mu<2$.

V.1.2 With broken time reversal symmetry

In the broken time reversal symmetry case, the topological invariant $\upsilon$, Eq. (8), is $\upsilon=-1$ for $-2w_{0}\cos\varphi_{0}<\mu<2w_{0}\cos\varphi_{0}$ and $\upsilon=1$ otherwise. For $\alpha>1.5$ it coincides with the presence or the absence of MZM, as shown in Fig. 9. For $1\leq\alpha\leq 1.5$ as in the previous case with time reversal symmetry, there is a coexistence of MZM and MEM together with regions where edge modes are absent. Finally in the regime with $\alpha<1$ the MZM seems to disappear and there is a wide region where MEM are present. However the larger the negative $\mu$ the greater is the $\alpha$ below which the edge modes disappear. We numerically discover a critical line, for $\alpha<1$, in the phase diagram, shown in Fig. 9, separating the region with massive edge modes from the region where there are no edge modes.

V.2 Long-range hopping and pairing, with and without TR symmetry

In this section we will consider the Kitaev chain with same long range in the hopping, $w_{\ell}=w_{0}d_{\ell}^{-\alpha}$ and in the pairing $\Delta_{\ell}=\Delta d_{\ell}^{-\alpha}$ terms, in the presence of time reversal symmetry ($\varphi_{\ell}=0$), or in its absence ($\varphi_{\ell}\neq 0$). Quite in general, the Hamiltonian can be written as it follows

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h_{0}(k)}=\sum_{\ell=1}^{\infty}w_{\ell}\sin\varphi_{\ell}\sin(k\ell), {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h_{z}(k)}=-{\mu}/{2}-\sum_{\ell=1}^{\infty}w_{\ell}\cos\varphi_{\ell}\cos(k\ell) and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h_{y}(k)}=\sum_{\ell=1}^{\infty}\Delta_{\ell}\sin(k\ell), which, taking the form $\varphi_{\ell}=\varphi_{0}+\varphi\,\ell$ for the phase, are

[TABLE]

where ${\textrm{Li}_{\alpha}(z)}$ is the polylogarithm of order $\alpha$ and argument $z$.

V.2.1 With time reversal symmetry

For the long-range case the expression for the winding number ${\sf w}$ defined in Eq. (7), because of the non-analycity of ${\cal H}(k)$, Eq. (22), at $k=0$ for $\alpha<1$, gives, as in the long-range pairing case, as a result ${\sf w}=\pm\frac{1}{2}$. However we can draw the critical lines as the $r\rightarrow\infty$ limit of Eq. (9), requiring $h_{z}(0)=h_{z}(\pi)=0$,

[TABLE]

defined by the Riemann zeta function and Dirichlet eta function. For $\alpha<1$, the winding number, from Eq. (7) and using Eqs. (23), (24), at $\varphi_{\ell}=0$, turns to be ${\sf w}=\frac{1}{2}$ for $\mu<\mu_{c1}$ and ${\sf w}=-\frac{1}{2}$ for $\mu>\mu_{c1}$. For $\alpha>1$ instead ${\sf w}=1$ for $\mu_{c1}<\mu<\mu_{c2}$ and ${\sf w}=0$ otherwise.

At this point a comment about the winding number is in order. As in the case of long-range pairing, the winding number for $\alpha<1$ (long range) seems to be half an integer number because of the divergence of ${\cal H}(k)$, in terms of polylogharithmics, at $k=0$, and in particular because the function $h_{y}(k)$ does not close varying $k\in[0,2\pi)$. Actually we observe that

[TABLE]

However if we perform the limit $k\rightarrow 0$ before $r\rightarrow\infty$, then $h_{y}(0)=h_{y}(\pi)=0$ therefore in this respect, by this regularization, the winding numbers come to be integer valued again, namely ${\sf w}=\pm\frac{1}{2}$ are replaced by ${\sf w}=1$, [math].

As far as the edge modes are concerned, for $\alpha\geq 1$ we get MZM for $\mu_{c1}<\mu<\mu_{c2}$ and no edge modes otherwise. From finite-size scaling analysis on the energy spectrum we get edge modes with zero masses (MZM) also for $\alpha$ slightly below $1$, for $-w_{0}\lesssim\mu\lesssim 2w_{0}\ln 2$. The regime with $\alpha\lesssim 0.6$ and $\mu<2w_{0}\eta(\alpha)$, instead, is characterized by massive edge modes (MEM), see Fig. 10.

V.2.2 With broken time reversal symmetry

We conclude this section considering the case of long-range hopping and pairing with broken time reversal symmetry, with $\varphi_{\ell}=\ell\varphi_{0}\neq 0$. By this choice for the breaking symmetry phase, the critical lines are

[TABLE]

so that the topological invariant is $\upsilon=-1$ for $\mu_{c1}<\mu<\mu_{c2}$ and $\upsilon=1$ otherwise. However we already experienced the fact that the topological invariants do not uniquely determine the presence of specific edge modes. By the analysis of the energy spectrum we get that for $\alpha>1$ the Majorana modes (MZM) are present in the topological phase denoted by $\upsilon=-1$. For $\alpha<1$ we can have MZM as well as MEM for $\mu<\mu_{c2}$, (see Fig. 11). Few examples of the finite-size scaling analyses are reported in Fig. 12.

VI Majorana wavefunctions

In this section we calculate explicitly the wavefunctions of the zero energy modes for the extended Kitaev model with generic $r$-neighbor hopping and pairing terms, both in the thermodynamic limit, by means of a generalized transfer matrix approach, and for a finite length chain, generalizing the discussion done for the first-neighbor model in Ref. Kitaev . For this purpose it is more convenient to adopt the Majorana mode representation. Let us introduce and write the Majorana operators in terms of the fermionic operators $a_{j}$ and $a_{j}^{\dagger}$

[TABLE]

which fulfill $\{c_{i},c_{j}\}=2\delta_{i,j}\;\forall\,i,j=1,\dots,L$ and $c_{j}^{\dagger}=c_{j}$. The Hamiltonian Eq. (II), therefore, can be rewritten in terms of the Majorana operators as follows

[TABLE]

We will now calculate the wavefunctions as solutions of the Bogoliubov-de Gennes equations (see Appendix A).

VI.1 Transfer matrix approach

We will consider for simplicity the case with time reversal symmetry ($\varphi_{\ell}=0$), so that Eq. (VI) reduces to

[TABLE]

Introducing the wavefunctions $\phi_{\epsilon_{n},j}$ and $\psi_{\epsilon_{n},j}$, related to the Bogoliubov coefficients which diagonalize $H$ ($\phi_{\epsilon_{n},j}=u_{n,j}+v_{n,j}$, $\psi_{\epsilon_{n},j}=u_{n,j}-v_{n,j}$), we get the following Bogoliubov equations (see Appendix A for details)

[TABLE]

For the most general case one has to solve Eqs. (81)-(84). Since we are interested in the zero energy states, we will consider the case where, for some ${n}$, we have $\epsilon_{n}=0$.

In this case the wavefunctions $\phi_{o,j}$ and $\psi_{o,j}$, related to zero energy, in Eqs. (VI.1), (35) decouple and we can solve the equations separately. Supposing $\Delta_{r}+w_{r}\neq 0$, Eqs. (VI.1), (35), for $\epsilon_{n}=0$, can be written in the following form

[TABLE]

after introducing the following transfer matrix

[TABLE]

where

[TABLE]

It is straightforward to notice that $t_{1}={\textrm{Tr}}(T)=\sum_{i=1}^{2r}\lambda_{i}$, is the trace of $T$, and $t_{2r}=-\textrm{Det}(T)=-\Pi_{i=1}^{2r}\lambda_{i}$, the determinant, where $\lambda_{i}$ are the eigenvalues of $T$.

By this approach one can write the wavefunctions evaluated at some point from its value at another point by applying several times the transfer matrix, namely applying $T$ to some power, say $j$, related to the space distance between the two points. What is relevant is therefore $T^{j}$, or, more conveniently, its diagonal form, so that $T^{j}$ can be written as $T^{j}=SD^{j}S^{-1}$ were $S$ diagonalizes $T$ ($D=S^{-1}TS$). The problem is reduced, therefore, to finding the eigenvalues ($D$) and the eigenstates ($S$) of $T$. In order to find the eigenvalues one has to write the charateristic polynomial $p_{2r}(\lambda)$ of the $2r\times 2r$ matrix in Eq. (38), and find the solutions of $p_{2r}(\lambda)=0$. One can easly prove that the polynomial $p_{2r}(\lambda)=\textrm{Det}(T-\lambda\mathbb{1})$ is such that

[TABLE]

and, therefore, by iteration, and making it equal to zero, one has to solve the following eigenvalue equation

[TABLE]

in order to find the $2r$ eigenvalues of $T$, $\lambda_{s}$ with $s=1,\dots,2r$. One can easily check that the corresponding eigenfunctions are $(\lambda_{s}^{2r},\lambda_{s}^{2r-1},\dots,\lambda_{s},1)^{t}$, which are the columns composing the matrix $S$. As a result, generic solutions of Eqs. (36), (37) can be written as the following linear combinations

[TABLE]

where the coefficients $c^{\phi}_{s}$ and $c^{\psi}_{s}$ are independent and $j=1,\dots,L$. A part of those coefficients can be fixed by imposing open boundary conditions. Let us call $n^{<}$ the number of $\lambda_{s}$ such that $|\lambda_{s}|<1$ and, analogously, $n^{>}$ the number of eigenvalues such that $|\lambda_{s}|>1$. Let us suppose the majority of $\lambda_{s}$ are smaller than $1$, namely if $n^{<}>n^{>}$, then the boundary conditions (open boundary condition and thermodynamic limit) require

[TABLE]

In this way we can have functions localized at the edges, $\phi_{o,j}$ on the left (small $j$) and $\psi_{o,j}$ on the right (large $j$). The eigenvalues greater than one in modulus will be discarted ($c^{\psi}_{s}=c^{\phi}_{s}=0$) to satisfy normalization conditions. The situation is reversed if, instead, $n^{<}<n^{>}$. In this case we have to impose the following boundary conditions

[TABLE]

verifying the existence of zero energy modes which will be localized at the edges, $\psi_{o,j}$ on the left and $\phi_{o,j}$ on the right.

In both the cases, the number of boundary conditions to impose for each wavefunction is $r$. The degrees of freedom to construct the wavefunctions in Eqs. (44), (45) are, therefore, equal to N=\big{(}\textrm{max}(n^{<},n^{>})-r\big{)}. This means that we can construct $N$ different linear combinations Eqs. (44), (45) at each edge, namely $N$ different zero-energy wavefunctions per edge. Since $\textrm{max}(n^{<},n^{>})\leq 2r$, we can get at most $N=r$ Majorana states per edge, consistently with the winding number analysis. A particular attention has to be paid for complex $\lambda_{s}$ since globally $\phi_{o,j}$ and $\psi_{o,j}$ have to be real valued by construction. This is also the reason of some oscillating behaviors of the wavefunctions, which come together with an overall exponential decay, since complex eigenvalues appear in pairs, therefore, $|\lambda_{s}|^{j}e^{i\theta j}+|\lambda_{s}|^{j}e^{-i\theta j}=2\,e^{j\ln|\lambda_{s}|}\cos(\theta j)$.

An example with $r=2$ is provided in Fig. 13, where two different $\phi_{o,j}$, two ortonormalized wavefunctions $\phi^{{}^{\prime}}_{o,j}$, $\phi^{{}^{\prime\prime}}_{o,j}$, both localized at the same edge, has been found for a set of parameters for which the winding number is ${\sf w}=2$. Analogously one can find two independent wavefunctions localized on the other edge of the chain.

Another simple example which shows the connection between this approach and the previous winding number analysis is the following. Let us consider for simplicity $r=2$ and $w_{1}=w_{2}=\Delta_{1}=\Delta_{2}=1$, or in terms of the parametrization in Eq. (2), $\alpha=\beta=0$ and $w_{0}=\Delta=1$. In this case Eq. (43) is simply $\lambda^{4}+\lambda^{3}+\frac{\mu}{2}\lambda^{2}=0$ whose solutions are all $|\lambda_{s}|<1$ for $0<\mu<2$, consistently with Eq. (17) for $r=2$.

Actually, the case $w_{\ell}=\pm\Delta_{\ell}$ has to be consider separately, as we will do in the next section, since in that case $T$ is singular and some of the eigenvalues are zero, requiring in some cases $\phi_{o,j}$ and $\psi_{o,j}$ to be Dirac delta functions, which are actually distributions, not really functions. The advantage in that case is that we can relax the condition of working in the thermodynamic limit, required for vanishing overlap between the wavefunctions localized at the two edges, which otherwise would hybridize spoiling the presence of unpaired Majorana modes.

VI.2 Special cases: Majorana modes in a finite-lenght chain

Let us now consider the model Eq. (VI.1) for the spacial case where $\Delta_{\ell}=w_{\ell}$ and, for simplicity $\mu=0$. In this case the Hamiltonian is simply given by

[TABLE]

The full model in Eq. (VI.1) can be sketched as in Fig. 14 where the case with $r=2$ as been shown.

Some of the couplings in Eq. (VI.1) disappear for $\mu=0$ and $\Delta_{\ell}=\pm w_{\ell}$.

For $\Delta_{\ell}=w_{\ell}$ the Hamiltonian Eq. (50) can be represented as in Fig. 15 where the simple case with $r=2$ is depicted.

In Eq. (50) two Majorana modes are missing, $c_{1}$ and $c_{2L}$, therefore we have one unpaired Majorana mode at each edge. This is consistent with winding number ${\sf w}=1$, Eq. (16) at $\mu=0$. Analogously, for $\Delta_{\ell}=-w_{\ell}$, the corresponding Hamiltonian would be H=-i\sum_{\ell=1}^{r}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{j=1}^{L-\ell}}w_{\ell}\,c_{2j-1}c_{2(j+\ell)}, in this case the unpaired Majorana operators, the ones which are missing in the Hamiltonian, are $c_{2}$ and $c_{2L-1}$.

Let us now consider the special case where $\Delta_{r}=w_{r}\neq 0$ and $\Delta_{\ell}=w_{\ell}=0,\;{\textrm{for}}\;\ell=1,\dots,r-1$. In this case the Hamiltonian Eq. (50) reduces to

[TABLE]

and we have $2r$ unpaired Majorana operators, $r$ at each edge, $c_{2j-1}$ with $j=1,\dots,r$ on the left side and $c_{2(L-j+1)}$ with $j=1,\dots,r$ on the right side. An example with $r=2$ is sketched in Fig. 16.

Defining the fermionic operators

[TABLE]

for $j=1,\dots,L-r$, Eq. (51) can be written as

[TABLE]

where $r$ highly non local fermions are missing, which can be defined as

[TABLE]

with $j=1,\dots,r$. Let us define the following parity operator

[TABLE]

Since these terms are missing in the Hamiltonian, then $[P,H]=0$, therefore $H$ and $P$ have the same eigenstates. Taking the states $|\pm\rangle$ such that $\tilde{a}_{j}|+\rangle_{j}=0$, $\tilde{a}_{j}|-\rangle_{j}=|+\rangle_{j}$ and $\tilde{a}_{j}^{\dagger}|+\rangle=|-\rangle_{j}$, defining $|\Psi_{0}\rangle=\left(\otimes_{j=1}^{r}|\pm\rangle_{L-r+j}\right)$, we have

[TABLE]

where $p$ is the number of single particle states $|-\rangle$ contained in $\left(\otimes_{j=1}^{r}|\pm\rangle_{L-r+j}\right)$ and $E_{0}=w_{r}(r-L)$ a constant which can be arbitrarily shifted to zero. Since we can construct $|\Psi_{0}\rangle$ choosing for each site $(L-r+j)$, $|+\rangle$ or $|-\rangle$, the degeneracy in energy is $2^{r}$, divided in two parity sectors, $2^{r-1}$ with odd parity and $2^{r-1}$ with even parity.

VII Conclusions

In this paper we studied the Kitaev chain, generalized by allowing for a longer range of the coupling terms. We considered the system when time reversal symmetry holds or is broken, characterizing the phases by topological invariants and by the spectrum properties, looking at the ground state energy and energy gaps. We show that when both the hopping and the pairing terms are extended to many neighbors, with time reversal symmetry, we can have many Majorana modes localized at the egdes. This finding is also confirmed by the explicit calculation of the edge modes by means of generalized transfer matrix approach useful to solve the proper set of Bogoliubov-de Gennes equations. The limit of strictly long-range couplings, with infinite neighbors, has also been considered, with and without time reversal symmetry, finding several phase diagrams, all characterized by the presence of both massless and massive edge modes.

Acknowledgements.

LD thanks SISSA for kind hospitality and acknowledges financial support from University of Padova through the project BIRD 2016. We thank N. Lo Gullo and L. Lepori for useful discussions.

Appendix A Bogoliubov transformation

The Hamiltonian in Eq. (II), in the fermionic representaion and real space, can be generally written in the matrix form

[TABLE]

and can be diagonalized by means of a unitary transformation

[TABLE]

getting the typical spectrum of a particle-hole symmetric Hamiltonian. On the other hand the Hamiltonian can be written in terms of Majorana operators as in Eq. (VI), which can be written quite in general as it follows

[TABLE]

where $A$ is a real antisymmetric matrix ($A_{ij}^{*}=A_{ij}=-A_{ji}$) such that

[TABLE]

where $\{\pm\epsilon_{n}\}$ are eigenvalues of ${\cal H}$ and $W$ is a real orthogonal matrix, $WW^{T}=W^{T}W=\mathbb{1}$ enciclopedia .

Under the actions of $U$ and $W$ on respectively fermionic and Majorana operators

[TABLE]

the Hamiltonian $H$ can be written in the following forms

[TABLE]

Assuming that $U$ and $W$ are canonical transformations we have $\left\{\tilde{a}_{n}^{\dagger},\tilde{a}_{m}\right\}=\delta_{n,m}$ and $\left\{\tilde{a}_{n},\tilde{a}_{m}\right\}=0,\,\forall\,n,m=1,\dots,L$, together with $\left\{b_{n},b_{m}\right\}=2\delta_{n,m},\,\forall\,n,m=1,\dots,2L$. The two sets of new operators are linked by

[TABLE]

The quadratic Bogoliubov-de Gennes Hamiltonian $H$, written in terms of fermionic operators, can be diagonalized by Bogoliubov transformation, which is represented by the matrix $U$. We can write $\tilde{a}_{n}$ as combination of the first set of operators $a_{j}$ by means of two sets of functions $u_{n,j}$ and $v_{n,j}$

[TABLE]

with $\sum_{j=1}^{L}\left(|u_{n,j}|^{2}+|v_{n,j}|^{2}\right)=1,\,\forall\,n=1,\dots,L$. Thus we can write $U^{\dagger}$ in the following way

[TABLE]

From Eq. (63) and Eqs. (66), (66), we can write

[TABLE]

or, viceversa,

[TABLE]

Now using the Heisenberg equations, derived from Eq. (64),

[TABLE]

expressing $b_{n}$ in terms of $c_{j}$ and using Eq. (61), we can write

[TABLE]

More specifically, using Eq. (VI) we get the following Bogoliubov equations of the wavefunctions for the extended Kitaev model

[TABLE]

If time reversal symmetry holds ($\varphi_{\ell}=0$) then we can choose $u^{*}_{n,j}=u_{n,j}$ and $v^{*}_{n,j}=v_{n,j}$, therefore $W_{2n,2j-1}=W_{2n-1,2j}=0$ while, calling $W_{2n-1,2j-1}$ the wavefunction $\phi_{\epsilon_{n},j}$ and $W_{2n,2j}$ the wavefunction $\psi_{\epsilon_{n},j}$

[TABLE]

Eqs. (81)-(84) reduce to Eqs. (VI.1), (35).