[Kwant] Plotting Andreev Spectrum of a 4 \pi Josephson Junction

Hi Kwant users ! I was trying to plot the 4 \pi periodic Andreev spectrum of a topological Josephson junction. I was using a finite TSC-QSH-TSC junction and trying to calculate the spectrum as suggested in ( https://mail.python.org/archives/list/kwant-discuss@python.org/message/BYVQE... ). I have attached my code and the plot I got from the simulation below. *Am I doing something wrong?* __________________________________________________________________________________________________ # InAs Sample # TSC-QSH-TSC Josephson Junction # 4 \pi periodicity of Andreev Spectrum import kwant import numpy as np import math ## import matplotlib matplotlib.use('Agg') ## # For plotting from matplotlib import pyplot from matplotlib import cm import scipy.sparse.linalg as sla # For matrix support import tinyarray tau_z = tinyarray.array([[1,0], [0,-1]]) tau_x = tinyarray.array([[0,1], [1,0]]) tau_y = tinyarray.array([[0,-1j], [1j,0]]) I4 = tinyarray.array([[1,0,0,0], [0,1,0,0],[0,0,1,0],[0,0,0,1]]) SzSz = tinyarray.array([[1,0,0,0], [0,-1,0,0],[0,0,-1,0],[0,0,0,1]]) SxSx = tinyarray.array([[0,0,0,1], [0,0,1,0],[0,1,0,0],[1,0,0,0]]) SySy = tinyarray.array([[0,0,0,-1], [0,0,1,0],[0,1,0,0],[-1,0,0,0]]) S2Sx = tinyarray.array([[0,0,0,1], [0,0,1,0],[0,1,0,0],[1,0,0,0]]) SySx = tinyarray.array([[0,0,0,-1j], [0,0,-1j,0],[0,1j,0,0],[1j,0,0,0]]) S3Sx = tinyarray.array([[0,0,0,-1j], [0,0,-1j,0],[0,1j,0,0],[1j,0,0,0]]) SzI = tinyarray.array([[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]]) ISz = tinyarray.array([[1,0,0,0], [0,-1,0,0],[0,0,1,0],[0,0,0,-1]]) SzSx = tinyarray.array([[0,1,0,0], [1,0,0,0],[0,0,0,-1],[0,0,-1,0]]) ISx = tinyarray.array([[0,1,0,0], [1,0,0,0],[0,0,0,1],[0,0,1,0]]) ISy = tinyarray.array([[0,-1j,0,0], [1j,0,0,0],[0,0,0,-1j],[0,0,1j,0]]) SyI = tinyarray.array([[0,0,-1j,0], [0,0,0,-1j],[1j,0,0,0],[0,1j,0,0]]) SxI = tinyarray.array([[0,0,1,0], [0,0,0,1],[1,0,0,0],[0,1,0,0]]) SySz = tinyarray.array([[0,0,-1j,0], [0,0,0,1j],[1j,0,0,0],[0,-1j,0,0]]) SzSy = tinyarray.array([[0,-1j,0,0], [1j,0,0,0],[0,0,0,1j],[0,0,-1j,0]]) SxSz = tinyarray.array([[0,0,1,0], [0,0,0,-1],[1,0,0,0],[0,-1,0,0]]) SxSy = tinyarray.array([[0,0,0,-1j], [0,0,1j,0],[0,-1j,0,0],[1j,0,0,0]]) # Define the tight-binding system def make_system(Bfields, W=100, L1 = 50 ; L2=100 ; L=150 ; Delta = 1.0): # Define InAs parameters A=37.0 ; B=-660.0 ; D= 0.0 ; M=-7.8 ; theta=Bfields*np.pi ; # superconducting phase bias syst = kwant.Builder() # Here, we are only working with square lattices a = 10 lat = kwant.lattice.square(a,norbs=8) # Define the scattering region #----------------------------------------------------------------------------------------------- #========================================================================================================= # Define left topological superconductor for j in range(0,W): for i in range(0, L1): # On-site Hamiltonian syst[lat(i, j)] = -4*D*np.kron(np.eye(4),tau_z)/a**2 - 4*B*np.kron(ISz,tau_z)/a**2 + M*np.kron(ISz,tau_z) + Delta*(np.cos(theta))*np.kron(np.eye(4), tau_x) + Delta*(np.sin(theta))*np.kron(np.eye(4), tau_y) if i > 0: syst[lat(i, j), lat(i - 1, j)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*np.kron(SzSx,tau_z)/(2*1j*a) if j > 0: syst[lat(i, j), lat(i, j - 1)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*1j*np.kron(ISy,tau_z)/(2*a) #--------------------------------------------------------------------------------------------------------------------- # Define Middle quantum spin-Hall system for i in range(L1, L2): syst[lat(i, j)] = -4*D*np.kron(np.eye(4),tau_z)/a**2 - 4*B*np.kron(ISz,tau_z)/a**2 + M*np.kron(ISz,tau_z) if i > 0: syst[lat(i, j), lat(i - 1, j)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*np.kron(SzSx,tau_z)/(2*1j*a) if j > 0: syst[lat(i, j), lat(i, j - 1)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*1j*np.kron(ISy,tau_z)/(2*a) #------------------------------------------------------------------------------------------------------- # Define right topological superconductor for i in range(L2, L): # On-site Hamiltonian syst[lat(i, j)] = -4*D*np.kron(np.eye(4),tau_z)/a**2 - 4*B*np.kron(ISz,tau_z)/a**2 + M*np.kron(ISz,tau_z) + Delta*np.kron(np.eye(4), tau_x) if i > 0: syst[lat(i, j), lat(i - 1, j)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*np.kron(SzSx,tau_z)/(2*1j*a) if j > 0: syst[lat(i, j), lat(i, j - 1)] = D*np.kron(np.eye(4),tau_z)/a**2 + B*np.kron(ISz,tau_z)/a**2 + A*1j*np.kron(ISy,tau_z)/(2*a) syst = syst.finalized() return syst #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B_min = 0.0 ; B_max = 8.0 ; npts=51 ; Bfields = np.linspace(B_min, B_max, npts) data=[]; for i in np.arange(len(Bfields)): syst = make_system(Bfields[i]) ham_mat = syst.hamiltonian_submatrix(args=[Bfields], sparse=True) # we only calculate the k lowest eigenvalues ev = sla.eigsh(ham_mat, k=4, which='SM', return_eigenvectors=False) data.append(ev) pyplot.figure() pyplot.plot(Bfields, data) pyplot.xlabel("Phase [in units of pi]") pyplot.ylabel("energy [t]") #pyplot.show() pyplot.savefig('InAs_Andreev_Spectrum_k=4.pdf') Thanks in advance.