Dear all,
i used a slightly modified "discretize.py" code of
https://kwant-project.org/doc/1/tutorial/discretize (see attached file)
and added the term
""" + kappa * kron(sigma_z,sigma_0) """
to the hamiltonian (splitting up the different spins). When checking the
unitarity of the scattering matrix at energy=0 by applying the function
"check_unitarity(matrix)" for
kappa=0 and kappa=10**(-9) the unitarity-error is increasing from order
10**(-13) to 10**(-6) !
Why is the error of $ S^\dagger S$ getting so high?
Thanks in advance!
Andreas Bereczuk
# Tutorial 2.9. Processing continuum Hamiltonians with discretize
# ===============================================================
#
# Physics background
# ------------------
# - tight-binding approximation of continuous Hamiltonians
#
# Kwant features highlighted
# --------------------------
# - kwant.continuum.discretize
import kwant
import kwant.continuum
import scipy.sparse.linalg
import scipy.linalg
import numpy as np
# For plotting
import matplotlib as mpl
from matplotlib import pyplot as plt
def qsh_system(a=30, L=60, W=3000):
hamiltonian = """
+ C * identity(4) + M * kron(sigma_0, sigma_z)
- B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
- D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ A * k_x * kron(sigma_z, sigma_x)
- A * k_y * kron(sigma_0, sigma_y)
+ kappa * kron(sigma_z,sigma_0)
"""
template = kwant.continuum.discretize(hamiltonian, grid=a)
def shape(site):
(x, y) = site.pos
return (0 <= y < W and 0 <= x < L)
def lead_shape(site):
(x, y) = site.pos
return (0 <= y < W)
syst = kwant.Builder()
syst.fill(template, shape, (0, 0))
lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
lead.fill(template, lead_shape, (0, 0) )
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
syst = syst.finalized()
#kwant.plot(syst,num_lead_cells=15,show=True)
return syst
def check_unitarity(matrix):
Mdagger=np.transpose(np.conjugate(matrix))
MdM=np.matmul(Mdagger,matrix)
absMdM=np.absolute(MdM)
shape=np.shape(MdM)[0]
identity=np.eye(shape)
errorMatrix=np.absolute(identity-absMdM)
maxError=errorMatrix.max()
error=10**(-7)
return np.allclose(absMdM,identity,atol=error),"max error= ", maxError
def SMatrix_properties(syst,energy,params):
SMatrix=kwant.smatrix(syst,energy,params=params)
S=SMatrix.data
print("S-Matrix unitarity= ",check_unitarity(S))
def analyze_qsh():
params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0 , kappa=0)
print("kappa= 0:")
syst = qsh_system()
kwant.plotter.bands(syst.leads[0], params=params,
momenta=np.linspace(-0.3, 0.3, 201), show=False)
plt.grid()
plt.xlim(-.3, 0.3)
plt.ylim(-0.05, 0.05)
plt.xlabel('momentum [1/A]')
plt.ylabel('energy [eV]')
#plt.show()
SMatrix_properties(syst,energy=0,params=params)
params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0 , kappa=10**(-9))
print("kappa= 10**(-9):")
syst = qsh_system()
SMatrix_properties(syst,energy=0,params=params)
def main():
analyze_qsh()
# Call the main function if the script gets executed (as opposed to imported).
# See <http://docs.python.org/library/__main__.html>.
if __name__ == '__main__':
main()
