[Kwant] Site ordering while plotting wavefunction

[Girish Sharma]

Hi,

I have attached a small code, trying to visualize Majorana wavefunctions at
the two end points of my 1D wire. The kwant.plotter doesn't really help, so
I'm trying basic pyplot. However, due to random ordering of sites, I do not
see the wavefunction peaking at the ends (see the last plot). Is there a
way out to do this ordering correctly in PYTHON?

Thank you for your time !

Best regards,
Girish

# Visualizing Majorana Fermions on a nanowire
# ================================================
#
# Physics background
# ------------------
#  s-wave proximity SC + Rashba SOC + magnetic field 

from matplotlib import pyplot
import kwant

# First, define the tight-binding system

sys = kwant.Builder()

# Here, we are only working with square lattices
a = 1
lat = kwant.lattice.square(a)
W = 1
L = 100

# Define the system hoppings 

t=1		# Hopping integral
mu=1		# chemical potential
h_x=1.5		# magnetic field in x	
Delta=0.5	# s-wave SC
alpha=2		# SOC strength

import tinyarray

# Define Pauli Spin Matrices
sigma_0 = tinyarray.array([[1, 0], [0, 1]])
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])

# Define Pauli Matrices for particle-hole
tau_0 = tinyarray.array([[1, 0], [0, 1]])
tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
tau_z = tinyarray.array([[1, 0], [0, -1]])

import numpy
kron=numpy.kron
import scipy.sparse.linalg as sla
	
for i in xrange(L):
    for j in xrange(W):
        sys[lat(i, j)]  = -mu * kron(sigma_0,tau_z) + h_x * kron(sigma_x,tau_0) + Delta * kron(sigma_0,tau_x)

	# Hopping in x-direction
        if i > 0:
              sys[lat(i, j), lat(i - 1, j)] = -t * kron(sigma_0,tau_z) -1j * alpha * kron(sigma_y,tau_z)

kwant.plot(sys)

# Finalize the system
sys = sys.finalized()

energies = []
    
# Obtain the Hamiltonian as a dense matrix
ham_mat = sys.hamiltonian_submatrix(sparse=True)


# we calculate the eigenvalues

evals,evecs = sla.eigsh(ham_mat, 10, which='SM',return_eigenvectors=True)

mylist = list(xrange(10))	
pyplot.figure()
pyplot.plot(mylist, evals,'.')
pyplot.xlabel("index")
pyplot.ylabel("energy")
pyplot.show()

# Visualize the wave functions in the system.

mylist = list(xrange(100))
prob=[]
pos=0; #visualize wavefunction of this eigenvalue index
prob= numpy.abs(evecs[::4, pos])**2 + numpy.abs(evecs[1::4,pos])**2 + numpy.abs(evecs[2::4, pos])**2 + numpy.abs(evecs[3::4, pos])**2

kwant.plotter.map(sys, prob,colorbar=False)
pyplot.figure()
pyplot.plot(mylist,prob)
pyplot.show()