Dear Sergey, I'm glad you found the solution, let me just add a remark about Kwant. Kwant excludes the modes corresponding to the eigenvalue lambda = 0, and projects all the degrees of freedom on the singular values of the hopping matrix. This leads to an improved stability, but makes it harder to extract the wave functions of the evanescent modes. Some times it isn't even possible to calculate the evanescent modes in a numerically stable way. We should really write down the details of the algorithm, but sadly we didn't get around doing that yet. Best, Anton On Tue, Aug 4, 2015 at 1:45 AM, Sergey <sereza@gmail.com> wrote:
Dear all, Let me post a reply to my own question about how to compute edge modes in kwant. The answer could be found in eq 3.15 of Michael Wimmer's thesis. I could not find a kwant ready-made implementation, so, I done it myself. One needs to compute the number of linearly dependent eigenvectors with eigenvalues \lambda = exp(i p z) inside the unit circle (i.e. decaying in a given direction). Below is a bit of my code that does the job.
def NEdgeModes(energy,kx,ky): H0=lead.cell_hamiltonian(args=[kx,ky]) H1=np.empty(H0.shape, dtype=complex) hop=lead.inter_cell_hopping(args=[kx,ky]) H1[:, : hop.shape[1]] = hop H1[:, hop.shape[1]:] = 0 Hsize=H0.shape[0] ##now construct A and B matrices for generalized eigenvalue problem A = np.empty((2*Hsize,2*Hsize), dtype=complex) B = np.empty((2*Hsize,2*Hsize), dtype=complex) A[:H0.shape[0],:H0.shape[1]]=np.zeros(H0.shape, dtype=complex) A[:H0.shape[0],H0.shape[1]:]=np.identity(Hsize, dtype=complex) A[H0.shape[0]:,:H0.shape[1]]=-H1.conjugate().transpose() A[H0.shape[0]:,H0.shape[1]:]=energy*np.identity(Hsize, dtype=complex)-H0 B[:H0.shape[0],:H0.shape[1]]=np.identity(Hsize, dtype=complex) B[:H0.shape[0],H0.shape[1]:]=np.zeros(H0.shape, dtype=complex) B[H0.shape[0]:,:H0.shape[1]]=np.zeros(H0.shape, dtype=complex) B[H0.shape[0]:,H0.shape[1]:]=H1 ev,vl=linalg.eig(a=A,b=B,right=True) u=vl.transpose()[:,:Hsize] evu=vl.transpose()[:,Hsize:] evmask=np.absolute(ev)<=1 ev0=ev[evmask] rank=np.linalg.matrix_rank(u[evmask,:], tol=1.0e-8) return ev0.shape[0]-rank
Best wishes, Sergey
On 31/07/15 17:11, Sergey wrote:
Dear Developers, I am a bit surprised that for my system with 16 sites and 9 hoppings between the unit cells I get lead.modes(E, args =[...] ) returning stabilized modes which are 9 vectors of length 9 . I understand that my hopping is degenerate and what I hoped to get from calling lead.modes(E, args =[...] ) is evanescent mode vectors of length 16 with corresponding complex momenta and some extra vectors spanning the non-propagating subspace.
Could you clarify how to get what I want, please.
Best wishes, Sergey