Analytical results for Green's functions of lattice fermions

TL;DR

This paper develops analytical methods using recurrence relations to compute Green's functions of lattice fermions, leading to new explicit formulas and the discovery of non-local corner states in topological superconductors.

Contribution

It introduces a recurrence relation approach for analytical Green's functions, providing explicit formulas and novel insights into topological phases.

Findings

Explicit Green's function for open Kitaev chain derived

Non-local fermionic corner states identified in 2D p-wave superconductor

Enhanced analytical tools for lattice fermion systems

Abstract

We present a further development of methods for analytical calculations of Green's functions of lattice fermions based on recurrence relations. Applying it to tight-binding systems and topological superconductors in different dimensions we obtain a number of new results. In particular we derive an explicit expression for arbitrary Green's function of an open Kitaev chain and discover non-local fermionic corner states in a 2D p-wave superconductor.

Equations

$$H_{\mathrm{tb}}=\sum _{i=1}^Nϵc_i^{†}{c}_i+\sum _{i=1}^{N-1}\gamma c_i^{†}{c}_{i+1}+{\gamma }^{\ast }{c}_{i+1}^{†}{c}_i,$$

$$g_{km}(\tau )=⟨T_{\tau }{c}_k(\tau )c_m^{†}(0)⟩,$$

$$c_i=\sqrt{\frac{2}{N+1}}\sum _{l=1}^N\sin (q_{l}i)d_l=S_{il}{d}_l,$$

$$q_l=\frac{\pi l}{N+1},1\le l\le N.$$

$$H=2\gamma \sum _{l=1}^N\cos (q_l)d_l^{†}{d}_l,$$

$$g_{km}(i{\omega }_n)=\sum _{p=1}^N{S}_{kp}\frac{1}{i{\omega }_n-E_p}{S}_{pm}.$$

$$g_{km}(i{\omega }_n)=\frac{2}{N+1}$$

$$g_{km}(i{\omega }_n)=\frac{U_{k-1}\left(\frac{i{\omega }_n}{2\gamma }\right)U_{N-m}\left(\frac{i{\omega }_n}{2\gamma }\right)}{\gamma U_N\left(\frac{i{\omega }_n}{2\gamma }\right)},$$

$$Z=\int \mathcal{D}[c^{†},c]e^{i{S}_1},S_1=\frac{1}{\beta }\sum _{i{\omega }_n}{c}^{†}(i{\omega }_n)\mathcal{A}c(i{\omega }_n),$$

$$c^{†}(i{\omega }_n)=(c_1^{†}(i{\omega }_n),\dots ,c_N^{†}(i{\omega }_n))$$

$$\displaystyle{\cal A}=\left(\begin{array}[]{cccccc}i\omega_{n}&-\gamma&0&\ddots&0&0\\ -\gamma&i\omega_{n}&-\gamma&\ddots&0&0\\ 0&-\gamma&i\omega_{n}&\ddots&0&0\\ \ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\ 0&0&0&\ddots&i\omega_{n}&-\gamma\\ 0&0&0&\ddots&-\gamma&i\omega_{n}\end{array}\right)\,.$$

$$g_{\mathrm{end}}(i{\omega }_n)=\underset{N\to \infty }{\lim }\frac{U_{N-1}\left(\frac{i{\omega }_n}{2}\right)}{U_N\left(\frac{i{\omega }_n}{2}\right)}$$

$$=\frac{i{\omega }_n}{2}+\sqrt{{\left(\frac{i{\omega }_n}{2}\right)}^2-1},$$

$$g_{11}^{-1}(N,i{\omega }_n)=g_0^{-1}(i{\omega }_n)-{\gamma }^2{g}_{11}(N-1,i{\omega }_n).$$

$$g_{\mathrm{bulk}}^{-1}(i{\omega }_n)=g_0^{-1}(i{\omega }_n)-2{\gamma }^2{g}_{\mathrm{end}}(i{\omega }_n).$$

$$H_{\mathrm{tb}}=\sum _{n=1}^N\sum _{m=1}^Mϵc_{n,m}^{†}{c}_{n,m}+\gamma \sum _{n=1}^{N-1}\sum _{m=1}^{M-1}{c}_{n,m}^{†}{c}_{n,m+1}$$

$$+c_{n,m+1}^{†}{c}_{n,m}+c_{n,m}^{†}{c}_{n+1,m}+c_{n+1,m}^{†}{c}_{n,m}.$$

$$Z=\int \left(\prod _{j=1}^M\mathcal{D}{c}_j^{†}\mathcal{D}{c}_j\right)e^{i{S}_2},$$

$$S_2=\sum _{i{\omega }_n}\sum _{j=1}^M{c}_j^{†}(i{\omega }_n)\mathcal{A}{c}_j(i{\omega }_n)$$

$$Z=\int \mathcal{D}{c}_M^{†}\mathcal{D}{c}_M{e}^{i{S}_2^{\prime }},$$

$$S_2^{\prime }=\sum _{i{\omega }_n}{c}_M^{†}(i{\omega }_n)({\mathbf{g}}_1^{-1}-{\mathbf{\Gamma }}^{†}{\mathbf{g}}_{M-1}\mathbf{\Gamma })c_M(i{\omega }_n).$$

$${\mathbf{g}}_M^{-1}={\mathbf{g}}_1^{-1}-{\mathbf{\Gamma }}^{†}{\mathbf{g}}_{M-1}\mathbf{\Gamma },$$

$$P_{M+1}$$

$$Q_{M+1}$$

$$\displaystyle\left(\begin{array}[]{c}P_{M}\\ Q_{M}\end{array}\right)=R^{M}\,\left(\begin{array}[]{c}P_{0}\\ Q_{0}\end{array}\right)\,,$$

$$\displaystyle R=\left(\begin{array}[]{cc}0&{\bf\Gamma}^{-1}\\ -{\bf\Gamma}^{\dagger}&{\bf g}_{1}^{-1}\,{\bf\Gamma}^{-1}\end{array}\right)\,.$$

$$\mathbf{\Gamma }={\mathbf{\Gamma }}^{-1}={\mathbf{\Gamma }}^{†}=1$$

$$\displaystyle R=\left(\begin{array}[]{cc}0&1\\ -1&{\cal A}\end{array}\right)\,.$$

$$\text{mbox}det\left[(1+{\lambda }^2)1-\lambda \mathcal{A}\right]={\lambda }^N{U}_N\left(\frac{1+{\lambda }^2}{2\lambda }-\frac{i{\omega }_n}{2}\right)=0.$$

$${\lambda }_{k1,2}=\left[\frac{i{\omega }_n}{2}-\cos \left(\frac{\pi k}{N+1}\right)\right]$$

$$\pm \sqrt{{\left[\frac{i{\omega }_n}{2}-\cos \left(\frac{\pi k}{N+1}\right)\right]}^2-1},$$

$$\displaystyle T=\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)$$

Full text

I Introduction

Fermionic lattice models are widely used not only as a purely theoretical tool but also as a basis for investigation and modelling of physical properties of real materials. Despite their relative formal simplicity – the Hamiltonians of many of them can be written down as bilinears of fermionic operators – many of them resist explicit analytical solution. The most prominent example is the Azbel-Hofstadter problem of lattice fermions subject to a magnetic field. Azbel64 ; Hofstadter76 Although it is possible to set up recurrence relations for relevant observables, an explicit analytical solution in a closed form has not yet been found. Remarkably, the presence of the gauge field complicates the progress making necessary the application of such advanced methods as Bethe ansatz.WIEGMANN1994495 ; Wiegmann94 ; ABANOV1998571 The situation is better for special constellations of fields, when some analytical solutions are possible, see e. g. [Hatsugai93PRB, ; Ueta1997, ]. But also in simpler field free situations analytical results for solutions of recurrence relations are rare and far between, vast majority of studies being concentrated on numerical treatment of the problem.LeeJoannopoulos1981 ; SurfaceToBulk2011 ; VelevButler_2004

If one restricts oneself to a system’s single particle Green’s functions (GFs) an explicit solution on the level of eigenstates is not always necessary. Here a formal functional integral for GFs can immediately be written down. However, its computation turns out to be very cumbersome. This kind of a functional integral is essentially a sum over all possible paths the particle can take during its evolution between two states. While in the absence of the gauge field the phase gathered along each individual path is simply related to the path length, in presence of the field it acquires a highly non-trivial dependence on the path geometry and topology. One attempt to take that into account is presented in [Belov1991, ]. However, the resulting expressions are complicated and very difficult to handle.

Recurrence relations for GFs were originally proposed in the series of works,Haydock1 ; Haydock_recursion and subsequently successfully used for numerical calculations in a great variety of set-ups, see e. g. [MacKinnon:1985rt, ]. In recent years this method became very popular in the field of topological insulators as it allows for fast and efficient band structure calculations in systems with any kinds of elementary cells and arbitrary coupling mechanisms between them, see e. g. [Chinese2006, ; KlichLong2014, ; DwivediChua16, ]. On the other hand, in many situations the recurrence relations can be solved analytically leading to compact and useful results for relevant physical quantities.DoornenbalStoof2015 ; MajoranaAK2016 The goal of our paper is twofold: first we apply this efficient technique to systems in gauge fields and second, we present new analytical results not only for systems with topologically non-trivial band structures but also for simple cubic lattices in different dimensions for different surface and bulk geometries.

The structure of our presentation is as follows: In Section II we explain our ideas on the simple example of a spinless 1D tight-binding chain. Among other things we present a new analytical result Eq. (8) for an arbitrary GF in terms of Chebyshev polynomials. Section III proceeds with 2D systems. First we set the stage and present a straightforward generalization of the 1D calculation, which yields an analog of Eq. (8) for the 2D case: the Eq. (43). After that we derive useful and simple expressions for the local density of states (DOS) at different positions in a lattice with open boundary. In Section III.2 we proceed with a system in a uniform magnetic field. Here our goal is to go beyond the seminal results of [Hatsugai93PRB, ] and [Belov1991, ] by deriving an explicit expression for the GF, which remains valid for arbitrary magnetic fields. For not too strong magnetic fields we derive analytical expressions for the GFs and assess their quality by comparison with the exact results. The subsequent two sections are devoted to topological superconductors. In Section IV we apply our method in order to derive Eq. (102), which is an explicit formula for an arbitrary GF of a Kitaev chain of finite length. This is one of the central results of the present work. In Section V we turn to 2D $p$-wave superconductor on a lattice of finite size. Using the previously developed approach we find non-local corner states akin to edge zero modes of open Kitaev chains and discuss their possible applications. Finally, in Section VI we present some previously unknown results for the 3D tight-binding lattice. A Conclusions section offers a short summary of our findings.

II Basics: 1D tight-binding chain

We begin our exposition with the simplest case – the spinless 1D tight-binding chain of length $N$ with the Hamiltonian

[TABLE]

where $\epsilon$ is the uniform energy on each site and $\gamma$ is the (in general complex) hopping amplitude between the sites. Our primary goal is the Matsubara GF of the form

[TABLE]

where $T_{\tau}$ is the imaginary time ordering operator. The most straightforward way to evaluate it is the direct diagonalization in terms of new operators $d_{l}$,

[TABLE]

where the momenta underlie the following quantization condition:

[TABLE]

The diagonalized Hamiltonian is then

[TABLE]

and the particle dispersion is obviously $E_{l}=2\gamma\cos(q_{l})$. In this representation the evaluation of the GF is straightforward and leads to

[TABLE]

Thus the result is

[TABLE]

The remaining sum can be evaluated and one obtains

[TABLE]

where $U_{k}(x)$ denote the Chebyshev polynomials of the second kind.gradshteyn2007 ; gogolin2013lectures ; PeresStauberDosSantos09 This expression holds for $k\leq m$, for $m>k$ we just have to interchange the indices. An alternative calculation can be performed using the functional integral formalism. The partition function is given by

[TABLE]

where

[TABLE]

are composite fields and the action kernel is given by an $N\times N$ matrix (from now on we concentrate on purely real $\gamma$, the case of generic tunnelling amplitude can be analyzed in exactly the same way)

[TABLE]

The matrix (8) of all possible GFs is just the inverse ${\cal A}^{-1}$. It can be calculated by the procedure proposed in [USMANI199459, ]. The resulting recurrence relations have solutions in terms of Chebyshev polynomials and one immediately obtains (8). We note in passing that similar methods can be used to treat systems with periodic boundaries. For basic results in 1D see Appendix A.

One practical application of the above result is the computation of the local density of states (DOS), which is found as the imaginary part of the local (at $k=m$) retarded GF. The latter is conveniently found from the analytically continued Matsubara GF via the substitution $i\omega_{n}\to\omega+i\delta$, where $\delta$ is a positive infinitesimal. Particularly interesting is the case of the edge site at $k=m=1$ or $k=m=N$ (from now on we use $\gamma$ as the energy unit),

[TABLE]

where the evaluation of the limit can be done in accordance with the procedure outlined in [Duran1999304, ] (for the proper analytical continuation see Appendix B). We note in passing that in different applications this kind of GF is also referred to as surface or boundary GF.

The same can be accomplished via the Dyson equation for the GF of the outmost chain site:

[TABLE]

Here $g^{-1}_{0}(i\omega_{n})=g^{-1}_{11}(1,i\omega_{n})=i\omega_{n}-\epsilon$ is the reciprocal of the Matsubara GF of an individual uncoupled chain site. In the limit $N\to\infty$ we can set $g_{\rm end}(i\omega_{n})=g_{11}(N,i\omega_{n})=g_{11}(N-1,i\omega_{n})$ and solve the corresponding equation. The GF of the bulk site (at $k=m=N/2$) in the limit of the infinitely long chain $N\to\infty$ can also be computed by a version of the above Dyson equation. In this situation $2g_{\rm end}(i\omega_{n})$ plays the role of the self-energy and replaces $g_{11}(N-1,i\omega_{n})$ in Eq. (19):111Obviously, the reversed procedure: from the bulk to the ‘surface’ GF is as simple, see also [VelevButler_2004, ].

[TABLE]

This equation can be considered to be the simplest version of the bulk-boundary correspondence, often considered especially in the context of systems with topologically non-trivial band structures.EssinGurarie

III 2D tight binding system

III.1 Zero field case

We construct a 2D system out of $M$ 1D systems of length $N$ arranged in parallel and coupled by the same matrix elements:

[TABLE]

The partition function for this system can be written down in terms of a functional integral over $M$ composite fields, which are this time arrays of objects (10):

[TABLE]

where

[TABLE]

Here ${\cal A}$ is as defined in (17) and $\boldsymbol{\Gamma}=-\mbox{diag}(\gamma,\dots,\gamma)$ is the diagonal matrix of rank $N$ coupling the chains. By a repeated integration over the fields $c_{1},\dots,c_{M-1}$ one obtains

[TABLE]

with the action

[TABLE]

Here ${\bf g}_{1}={\cal A}^{-1}$ denotes the matrix of all GFs for an individual 1D chain with length $N$, given in Eq. (8). Thus the GF matrix for the sites at the edge of the system is found from the recurrence relation

[TABLE]

with the initial value ${\bf g}_{0}=0$. This is a direct generalization of the relation (19) and also has the form of a Dyson equation.

As is shown in Ref. [MajoranaAK2016, ] Eq. (21) allows for an explicit solution in terms of matrix polynomials of Chebyshev type.FaberTichy2010 Here we take another route and write the unknown GF as a quotient of two matrices ${\bf g}_{M}=P_{M}\,Q_{M}^{-1}$. Then the recurrence relation can be split into two:

[TABLE]

which can immediately be solved in terms of powers of a matrix $R$:

[TABLE]

where

[TABLE]

Obviously the initial conditions are $P_{0}=0$ and $Q_{1}=1$. Let us now assume that $\gamma$ is purely real (this assumption considerably simplifies calculations and is not restrictive in any way) and our energy unit. Then

[TABLE]

and

[TABLE]

We would like to find $T$ and a diagonal $R_{0}$, so that $R=T\,R_{0}\,T^{-1}$. To that end we need the eigenvalues and -vectors of the matrix $R$. Its characteristic equation for the eigenvalues $\lambda$ reads

[TABLE]

This equation can be solved using the trigonometric representation of the Chebyshev polynomials. The $2N$ different solutions are given by

[TABLE]

where $1\leq k\leq N$. Therefore $R_{0}=$diag$(\lambda_{j1},\lambda_{j2})$. The matrix $T$ can be written down in terms of eigenvectors $v_{1,\dots,N}$ of the action matrix ${\cal A}$ for a 1D chain:

[TABLE]

with

[TABLE]

$v_{j}$ constitute a self-inverse symmetric matrix $V$, which is essentially a square root of a unity matrix,

[TABLE]

In fact, this result is already contained in the diagonalization transformation (3). Computation of powers of $R$ is now straightforward and is just

[TABLE]

where

[TABLE]

Due to the special property of the recurrence initial conditions only the right column of the result (39) is important for the GF computation. Plugging this back into (27) and computing the ratio of $P_{M}\,Q_{M}^{-1}$ we obtain the final result:

[TABLE]

In Eq. (41) one can immediately recognize the 1D result (7). On the other hand, from (7) and (8) follows the identity

[TABLE]

where

[TABLE]

Therefore the GF between the sites $(i,p)$ and $(j,q)$ is

[TABLE]

As an application we compute the DOS at different points in the lattice, see Fig. 1: (a) at the corner site of the lattice; (b) at the site on the edge of the system far away from the corners, which we call edge bulk (eb) site; (c) at the bulk far away from the edges. While an analytical result for (c) exists and is reported in [2D_DOS_old, ], the situations (a) and (b) have not yet been considered.

First we concentrate on the corner site $(1,1)$ of our lattice for $M\to\infty$ and $N\to\infty$. Taking the latter limit amounts to a replacement of the quotient of the polynomials in (41) by (II) with appropriate arguments. Using the result of the analytical continuation given in Appendix B and replacing the sum in (41) by an integral we then obtain the following result for positive energies $\omega>0$:

[TABLE]

With the help of a similar procedure we can compute the DOS at a site $(1,N/2)$ in the middle of the system edge, at the edge bulk site. Here in the limit $N\to\infty$ we obtain

[TABLE]

Finally, the genuine bulk GF can be computed using the version of the bulk-boundary correspondence condition following from Eq. (21):

[TABLE]

This equation describes a 1D system of length $N$ with the GF ${\bf g}_{1}$ (the chain with open circles on Fig. 1), which is coupled to two identical systems of sizes $N\times M$ with edge GFs ${\bf g}$.222These systems can, of course, in general be different. Since

[TABLE]

where $C=\mbox{diag}(i\omega_{n}-2\cos[\pi k/(N+1)])$ we thus obtain

[TABLE]

Using this result we derive an alternative expression for the edge bulk DOS of an infinitely large system:

[TABLE]

which yields exactly the same result as Eq. (45). For the genuine bulk DOS we then obtain the known result 2D_DOS_old

[TABLE]

The expressions (44), (45), (46), (47) are valid for positive energies $0<\omega<2$, for $\omega>2$ (outside of the band) the DOS is zero in all three cases. We would like to remark that all four integrals can be expressed in terms of elliptic functions. We abstain from doing so as it does not produce any added value.

III.2 2D lattice in magnetic field

In order to include magnetic field into the model (III.1) we use Landau gauge. We apply the Peierls substitution in the form depicted in Fig. 2.Peierls33

In this kind of geometry 1D systems of length $M$ are coupled by ‘bare’ tunnelling amplitudes $\gamma$, which do not contain magnetic field dependent phases. On the other hand, the tunnelling amplitude within each of such 1D chains is dressed by factors $e^{\pm im\phi}$, where $0<m<M-1$ is the chain index. The phase is defined as $\phi=Ba_{0}^{2}/\Phi_{0}$, where $B$ is the field magnitude, $a_{0}$ is the lattice constant and $\Phi_{0}=h/2e$ is the magnetic flux quantum.

The GF for particles in each of the $M$ chains is in analogy to (17) given by the inverse of the corresponding $N\times N$ action matrix:

[TABLE]

Unsurprisingly, the eigenvalues $\epsilon_{k}=i\omega_{n}/2-\cos[\pi k/(M+1)]$ ($1<k<M$) of this matrix do not depend on $\phi$ and the matrix of its eigenvectors is given by a still self-inverse matrix

[TABLE]

so that ${\cal A}_{m}=U\,B\,U$, where $B=\mbox{diag}(2\epsilon_{n})$. For the GF, which is a direct generalization of (8) one then obtains

[TABLE]

Such 1D systems are coupled to each other by the same matrices $\boldsymbol{\Gamma}$ as in the field free case. Therefore the recurrence relation for the GFs of the edge chain reads

[TABLE]

Its solution can be constructed in terms of a matrix $R_{s}$ defined as

[TABLE]

In this notation the above recurrence relation reads

[TABLE]

where as in the previous subsection ${\bf g}_{M}(i\omega_{n})=P_{M}Q_{M}^{-1}$, $P_{0}=0$ and $Q_{0}=1$. Its solution obviously is

[TABLE]

We use the following substitution:

[TABLE]

where $R_{0}$ is the matrix (33) of the field-free case and

[TABLE]

With this notation the above matrix product reads

[TABLE]

This representation is very useful for expansion in small fields as $D_{s+1}$ and thus $\epsilon_{s+1}$ are objects of the order $s\phi$. Therefore the expansion to the order $(M\phi)^{2}$ is given by

[TABLE]

In the next step we keep only terms linear in $\epsilon_{s}$. Since

[TABLE]

where ${\bf U}_{k}=\mbox{diag}\,U_{k}(\epsilon_{j})$ denotes a diagonal matrix containing Chebyshev polynomials of the argument $\epsilon_{j}$, $1<j<N$ as defined in Eq. (40), using the identity $U_{M-1}U_{M-s}-U_{M-s-1}U_{M}=U_{s-1}$ we obtain the following result for the GF on the edge of our system:

[TABLE]

The form of this correction allows for an interesting and useful interpretation. In accordance with Eqs. (8) and (43) one can consider the factors $V\,{\bf U}_{s-1}\,{\bf U}^{-1}_{M}\,V$ as being the GFs for the lattice nodes located in the rows $m=s$ and $m=M$, see Fig. 2. So the first factor from the right describes particle propagation from the row $m=M$ to the row $m=s$ where the particle ‘feels’ the magnetic field, picks up the factor $D_{s}$ and after that propagates back to the row with index $M$ thereby correcting the field free result ${\bf g}(i\omega_{n},0)$. According to this scheme the second order contribution can be found to be given by the following expression:

[TABLE]

Here the particle travels to the row $s^{\prime}$, picks up the phase induced by $D_{s^{\prime}}$, travels to the row $s>s^{\prime}$, picks up the second phase due to $D_{s}$ and after that returns back. In the similar way one can construct corrections of arbitrary order.

III.3 Edge states: the local density of states

It is known that in 2D systems subject to strong magnetic fields there are gapless edge states. It is interesting to recover them using the just developed technique. The simplest quantity is the local density of states (DOS), the spacial dependence of which is exemplarily plotted in Fig. 3.

A direct comparison of the edge DOS in the 2D case in presence of magnetic field with the one of a 1D system indicates the presence of edge channels with approximately 1D geometry.

III.4 Edge states: the edge currents

As we have seen in the previous subsection the local DOS is not sensitive to chirality of the edge states. In order to access this information a direct computation of particle currents is more appropriate. One can for example consider the current flowing between the sites $(n,M)$ and $(n+1,M)$ at the right outmost edge of a sample with dimensions $N$ and $M$. The corresponding operator is given by

[TABLE]

This quantity can be computed in the following way. Let us consider two 2D systems in the same magnetic field. The operators of one of them are $c_{n,m_{1}}$ and the operators of the other one $d_{n,m_{2}}$, where $1\leq n\leq N$ and $1\leq m_{1,2}\leq M_{1,2}$. We gauge the field in such a way that the (for definiteness lower) edge $m_{1}=0$ of the first system does not carry phase factors and they grow in positive direction for growing index $m_{1}$ up to the value $(M_{1}-1)\phi$ in the opposite edge of the system. On the other hand, the (upper) edge of the second system carries the phase factor $-\phi$ and it grows into negative values up to $-M_{2}\phi$ on the opposite edge of this subsystem.

Let ${\bf g}$ be the GF of the sites on the lower edge of the system $1$ and ${\bf h}$ be the GF of the upper edge of the second system. Now we couple the systems by tunnelling. Then the effective action has the form

[TABLE]

where the notation $({\bf c}^{\dagger},{\bf d}^{\dagger})$ stands for the composite field $(c_{1}^{\dagger},c_{2}^{\dagger}\dots c_{N}^{\dagger},d_{1}^{\dagger},d_{2}^{\dagger}\dots d_{N}^{\dagger})$ and $\Gamma$ is a unit matrix times $-\gamma$. We are interested in the GF between the neighboring lattice sites and we choose them to be of the kind $\langle T_{\tau}\,c_{n}(\tau)\,d^{\dagger}_{n}(0)\rangle$ and $\langle T_{\tau}\,d_{n}(\tau)\,c^{\dagger}_{n}(0)\rangle$. The difference between the two is precisely the current according to the definition (77). These expectation values can be found via matrix inversion of the above action. Taking the difference of the off-diagonal components of the inverse matrix we obtain an array of currents

[TABLE]

This expression is odd with respect to exchange ${\bf g}\leftrightarrow{\bf h}$. This automatically yields currents of opposite signs through the links located at the same distance from the middle symmetry axis of the sample. The opposite signs for the currents on the opposite sample edges follow. Moreover, in the absence of the magnetic field the matrices ${\bf g}$ and ${\bf h}$ commute and the net current through the links vanishes. This is due to the fact that both subsystems are diagonalized by the same transformation.

Numerical evaluation of the above expression is not difficult and one can conveniently discuss all features of the currents in the sample. Among other things one immediately verifies that the currents decay exponentially with the distance from the sample edge, see Fig. 4. We also find that there is no noticeable net field-driven depletion of charge even in the case of very strong fields.

It turns out that the edge currents computation can be very conveniently performed using the perturbative expansion presented in Section III.2. To the lowest order we can write ${\bf g}\approx{\bf g}_{0}+{\bf g}_{1}$ and ${\bf h}\approx{\bf h}_{0}+{\bf h}_{1}$ and use the result (III.2). Defining ${\bf f}:=\mathbb{1}-\gamma^{2}{\bf g}_{0}{\bf h}_{0}$ and taking advantage of the identity $\left[{\bf g}_{0},{\bf h}_{0}\right]=0$ one obtains

[TABLE]

where we defined $J_{i\omega_{n}}$ as the energy-resolved current with the property $J=\frac{1}{\beta}\sum_{i\omega_{n}}J_{i\omega_{n}}$, as well as $G_{0,1}=V{\bf g}_{0,1}V$, $H_{0,1}=V{\bf h}_{0,1}V$ and $F=V{\bf f}V$.

Assuming $\gamma$ to be the energy unit as before one arrives at an explicit expression for the current to first order in $e^{iM_{1}\phi}-1$ or $e^{iM_{2}\phi}-1$, whichever is larger:

[TABLE]

where $\epsilon_{k}$ were defined in (40) and the very last term denotes the term identical to the first one up to the indicated substitutions. In the last computation step we have performed an index shift of $s$, which does not affect the index of the Chebyshev polynomials.

Even in the expansion of this order one can see many features imposed on the system by the magnetic field. A comparison between the full current and the approximation can be found in Fig.4 and Fig. 5. As expected, the current decays exponentially with the distance from the sample edge. It turns out, that the perturbative approach works surprisingly well yielding a good approximation for the edge current up to the depths of about 10% of the sample size, see Fig. 4.

IV Kitaev chain model

Recently Kitaev chain model moved in the focal point of research as it represents one of the simplest realizations of non-local Majorana edge states.Kitaev2001 ; QuantumComputingReview2008 ; Alicea2012 Its Hamiltonian is given by

[TABLE]

where $H_{\text{tb}}$ is defined in (1). $\Delta$ is the gap parameter and $\phi$ is the superconducting phase. The action can be written in the form (9) after the introduction of composite fields of the form $c^{\dagger}(i\omega_{n})=(c_{1}^{\dagger}(i\omega_{n}),c_{1}(i\omega_{n}),\dots,c_{N}^{\dagger}(i\omega_{n}),c_{N}(i\omega_{n}))$. The action kernel is

[TABLE]

where the inter-site coupling matrix is

[TABLE]

and $D_{0}^{-1}=\mbox{diag}(i\omega_{n}-\epsilon,i\omega_{n}+\epsilon)$ is the action of a single detached site. The diagonal $2\times 2$ Nambu blocks of the inverse of ${\cal A}_{N}$ yield the local GFs of a respective site, on its diagonal are the electron and hole GFs, while the anomalous GFs are on the off-diagonals:

[TABLE]

where we define the local Matsubara GFs in Nambu representation as

[TABLE]

The non-local GFs are defined in the same way with the pair of $k$ indices replaced by the indices of the sites between which the respective GF is considered. Let us denote by $\widetilde{D}_{N}$ the GF of the kind (96) at the left outermost site of a chain with length $M$. Then the following recurrence relation holds:

[TABLE]

It can, of course, be recast into the Dyson equation as

[TABLE]

After making the substitution

[TABLE]

one then obtains a matrix-valued three-point recurrence relation for a new variable $\widetilde{P}_{N}$,

[TABLE]

It can be considered to be a matrix generalization of the Chebyshev polynomials of the second kind. Using a similar approach we can compute the GF on the outmost right site of the chain. The corresponding identities are obtained by the exchange $\Lambda\leftrightarrow\Lambda^{\dagger}$. We shall denote these GFs by $D_{N}$ and the respective Chebyshev polynomials by $P_{n}$.

Now we compute the local GF at the site $1<k<N$ of the chain. Obviously, it cuts the chain into two pieces: (i) with length $k-1$ to the left and (ii) with length $N-k$ to the right of the site $k$. That is why $D_{k-1}$ and $\widetilde{D}_{N-k}$ are the corresponding self-energies and

[TABLE]

Using similar procedures one can derive an explicit expression for the GF between arbitrary sites $k$ and $m$. It is given by

[TABLE]

For the detailed derivation see Appendix C. One special case: $k=1$, $m=N$ is particularly interesting as the respective GF is responsible for the transport properties of the chain. Here the expression is very appealing and concise:

[TABLE]

It is not difficult to show that in the scalar case the last two results immediately reproduce the corresponding formula for the tight-binding chain (8). One of the applications for that is the derivation of the effective action for end Majoranas used in Refs. [GolubHorovitz2011, ; Weithofer2014, ].

The most important advantage of these results is that the computation of the Chebyshev matrix polynomials itself only requires matrix multiplications. Only the very last steps in (102) and (103) require matrix inversions. This is different from directly using the recurrence relation (104), which, being a generalization of a continuous fraction to matrices, requires a matrix inversion in each step. Needless to say, it is also more efficient than the direct matrix inversion of (90).

Just as in the case of the simple tight binding chain one can produce a bulk-boundary correspondence relation using the Dyson equation (98).

[TABLE]

with the only difference that in the present situation there are two different edges coupled to a bulk site: a ‘right’ and a ‘left’ one.

In order to access the end site GF we can alternatively follow the route discussed in the previous section and define new matrices $P_{N}$ and $Q_{N}$, similar to those in Eq. (22), so that $D_{N}=P_{N}Q_{N}^{-1}$. Then the recurrence relation (97) is solved by

[TABLE]

where $P_{1}=D_{0}$ and $Q_{1}=1$. Further progress is made by diagonalizing the matrix $R$. Its eigenvalues are

[TABLE]

Eigenvectors of $R$ can be written down as $({\bf u},{\bf v})^{T}$. Obviously,

[TABLE]

From the requirement

[TABLE]

we obtain

[TABLE]

Setting $v_{1}=1$ we obtain a set of 4 different vectors ${\bf v}(\lambda_{1,2,3,4})$, from which we compute ${\bf u}(\lambda_{1,2,3,4})$ using (113). Then

[TABLE]

is the matrix which diagonalizes $R$. Thus we obtain

[TABLE]

In this way we obtain an analytical solution for the local GF in a Kitaev chain of finite length $N$.

V 2D $p$-wave superconductor

Kitaev chain model dealt with in the previous section is a 1D version of the more general $p$-wave superconductor models. In 2D it can be understood as a stack of Kitaev chains coupled by superconductor pairing (see e. g. [bernevig2013topological, ]):

[TABLE]

where $H_{\rm tb}$ is defined in Eq. (III.1). The system can be considered to be build up from $M$ Kitaev chains of length $N$, which are coupled by matrices $\boldsymbol{\Gamma}$:

[TABLE]

Then the recurrence relation for the GF of the edge row is formally equivalent to that given in Eq. (21):

[TABLE]

where ${\bf G}$ is the GF of an individual uncoupled Kitaev chain computed in (102). The solution of this recurrence relation can again be performed using yet another set of Chebyshev matrix polynomials of the second kind.

As an application we compute the energy spectra of systems of different size. As can be seen from

Fig. 6 there is a pronounced even-odd effect. While for even $N=M$ there are no zero modes, for odd $N=M$ there is always a double-degenerate energy level at $E=0$. For more generic lattice sizes the zero modes exist whenever both $N$ and $M$ are odd. The case of $M=1$ corresponds to an ordinary Kitaev chain, in which the zero modes lie at precisely $E=0$ for odd $N$ and approach zero energy with growing even $N$. A similar phenomenon takes place in the present case – also in the case of either of $N$ or $M$ being even (or both) the eigenenergies tend towards $E=0$ with growing lattice sizes. There is, however, a fundamental difference between these lattices and those with both $N,M$ being odd. While the former are fourfold degenerate, the latter are always doubly degenerate.

The present 2D case is obviously different from the edge state in the Kitaev chain. It is instructive to investigate the spatial distribution of the DOS in order to find out whether the zero modes are localized. As we have shown above in Fig. 1 there are three distinct points of the lattice: (a) corner site, (b) edge bulk site and (c) the true bulk site. In Fig. 7 we plot the corresponding DOS. Surprisingly, the zero energy states are indeed localized at the lattice corners and the respective DOS falls off exponentially with the distance from the corner. We would like to stress that for large $N$ and $M$ there are always states in vicinity of $E=0$. But only the odd-odd constellation possesses a true zero mode, which is doubly degenerate. In all other cases there are multiple levels (at least four) in the vicinity of $E=0$.

In the odd-odd case one is confronted with a fermionic state, which is highly delocalized between the four corners of the lattice. This is very similar to the end states in the open Kitaev chain.Kitaev2001 However, in the 2D case the corner states are not Majorana fermions. Nonetheless, they are as well perfectly suitable to be used as qubit states in the context of quantum information technology.

VI 3D tight binding lattice

The recurrence relation method can be straightforwardly generalized to lattices of higher dimensions. Here we consider a 3D cubic tight binding sample with dimensions $M\times N\times K$. We recover all known results and generate a number of new ones. The expression for the single-particle GF between the sites with coordinates $(jpq)$ and $(j^{\prime}p^{\prime}q^{\prime})$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

For the DOS at the corner site of the lattice we then obtain the following result

[TABLE]

whereby the integration domain is fixed by the requirements $|y|<1$, $|z|<1$ and $|\omega/2-y-z|<1$. The edge bulk site is the one on the edge of the sample far away from the corners

[TABLE]

with the same integration domain. In the bulk of the face far away from the edges we obtain

[TABLE]

And, finally, in the bulk of the lattice one finds

[TABLE]

To the best of our knowledge Eqs. (122) and (123) represent new results and are complementary to e. g. those of [Guttmann2010, ]. All remaining integrals can be rewritten in terms of elliptic integrals. We refrain from that though as it does not produce any added value.

VII Conclusions

We revisit the recurrence relation method (also referred to as transfer matrix method) for band structure calculations of lattice models and apply it for the computation of Green’s functions (GF). We show a number of analytical solutions for conventional lattices in different dimensions. While for a simple 1D tight-binding system every single GF can be written down as a rational function of Chebyshev polynomials of the second kind, in higher dimensions or for systems with such non-trivial structure as a $p$-wave superconductor or ones with spin-orbit coupling the resulting expressions for the GFs are given by functions of matrix-valued Chebyshev polynomials. Using this concept we derive an explicit expression for any kind of GF for a Kitaev chain, see Eq. (102). Even though our results require numerical calculations, their efficiency is vastly superior to all existing methods as they require only small number of matrix inversions, whatever the system size.

We generalize the method for lattices in external fields and for corresponding GFs derive analytical formulas in closed form. By an explicit computation of the density of states and energy-resolved particle currents we show how in a 2D tight binding lattice subject to a magnetic field a 1D chiral edge state is formed. Furthermore, we develop a perturbative approach in order to construct simpler analytical solutions, which adequately describe the properties of edge states for not too strong fields.

Application of our method to a 2D $p$-wave superconductor on a lattice reveals its very interesting energy level structure. It turns out, that for lattices with odd length and width there is always a doubly degenerate zero energy state, which is non-local and spread between the four corners of the lattice. That is confirmed by an explicit calculation of the spatial dependence of the respective density of states. This phenomenon is very similar to end states in open Kitaev chains. However, the emergent fermionic state is not of Majorana type. Nonetheless, its non-locality is very advantageous for future quantum information technology applications, not least due to its explicit higher dimensionality, which entails a better experimental feasibility.

In addition we derive numerous analytical results for densities of states of conventional tight-binding lattices in different dimensions and different spatial locations in lattices with open boundaries.

Acknowledgements.

AK is supported by the Heisenberg Programme of the Deutsche Forschungsgemeinschaft (Germany) under Grant No. KO 2235/5-1.

Appendix A

With the periodic boundary condition the action matrix is slightly modified:

[TABLE]

Its determinant and the inverse is given by

[TABLE]

On the other hand the matrix

[TABLE]

straightforwardly diagonalizes ${\cal A}$. That is why the GF is also given by the following expression:

[TABLE]

where $\epsilon_{k}=-2\gamma\cos(2\pi k/N)$. From the equality $g_{km}(i\omega_{n})=({\cal A}^{-1})_{km}$ then follows an interesting result for the trigonometric sum in (133).

Appendix B

The retarded GF obtained from the Matsubara GF in (II) can be found to be

[TABLE]

As a result the DOS is only nonzero within the band $|\omega|<2\gamma$ and has the expected half-ellipsoidal form

[TABLE]

In the bulk DOS one immediately recognizes the van Hove singularities,

[TABLE]

Appendix C

We first outline the calculation of the end-to-end GF. Here we use the chain contraction procedure proposed in MajoranaAK2016 . One starts with a partition function generated by the action (90) and subsequently integrates out all fermionic fields up to those describing the end sites. We start with

[TABLE]

Then, after integrating out the fermions of the second site (counted from the left) the action kernel is given by a matrix of reduced dimensions:

[TABLE]

Continuing this lattice contraction we finally obtain the $2\times 2$ action

[TABLE]

where for $F_{n}$, $T_{n}$ and $\widetilde{T}_{n}$ we have the following recurrence relations:

[TABLE]

with the initial conditions $F_{1}=D_{0}$, $T_{1}=\Lambda$ and $\widetilde{T}_{1}=\Lambda^{\dagger}$. So $F_{n}$ satisfy the recurrence relation for the local GF on the outmost right site $D_{n}$, see Section IV. Now we use the formula for the block matrix inversion

[TABLE]

where $A,B,C,D$ are arbitrary non-singular quadratic matrices. Applying this to invert (163) we immediately recognize, that

[TABLE]

i. e. it is equal to the local GF of the outmost left site of the chain. That is why for the object $G_{1N}$ we obtain the relation

[TABLE]

In order to access the $k$-to-$m$ GF for $m>k$ we follow the following strategy. First one integrates out all sites with indices $j<k$ and $j>m$. This process can be understood as outer chain contractions. Then one is confronted with the effective action

[TABLE]

Our task is now the calculation of the end-to-end GF for this action. To this end we employ the approach used above, which we call inner contractions. As a result, for the end-to-end effective action kernel we get

[TABLE]

Inversion of this expression yields for the off-diagonal component the value

[TABLE]

We would like to remark that this kind of procedure can be used for inversion (and, of course, calculation of the determinants) of block-tridiagonal matrices even in the general case of arbitrary $\Lambda$, $\Lambda^{\dagger}$ and $D_{0}$ matrices. Our approach is different from that presented in e. g. [Molinari97, ].