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.
Hi Subhadeep, At a glance, everything looks correct, except the eigenvalues aren't sorted. Please check out this blog post that explains the problem and possible solutions: https://quantumtinkerer.tudelft.nl/blog/connecting-the-dots/ Best, Anton On Thu, 14 Dec 2023 at 11:46, Subhadeep Chakraborty <sc20rs001@iiserkol.ac.in> wrote:
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.
Thanks for the suggestion. I'll look into this. Reagrds, Subhadeep On Thu, Dec 14, 2023 at 4:43 PM Anton Akhmerov <anton.akhmerov+kd@gmail.com> wrote:
Hi Subhadeep,
At a glance, everything looks correct, except the eigenvalues aren't sorted. Please check out this blog post that explains the problem and possible solutions: https://quantumtinkerer.tudelft.nl/blog/connecting-the-dots/
Best, Anton
On Thu, 14 Dec 2023 at 11:46, Subhadeep Chakraborty <sc20rs001@iiserkol.ac.in> wrote:
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.
participants (2)
-
Anton Akhmerov -
Subhadeep Chakraborty