Hello Kwant,
The problem I am trying to solve relies on being able to calculate the
eigenphases of the scattering matrix for different energies, more precisely
it requires the elements of the scattering matrix to be smooth as a
function of energy.
By examining the elements of the scattering matrix as a function of energy,
I have noticed that the elements of the scattering matrix are not smooth
and there exists a sign problem.
The code I have included below is for a highly(!!) simplified example,
where the scattering region is simply two leads attached to a single
square. I took the base code from the kwant tutorial page to reduce any
possible errors. I have taken a small range of energies and plotted the
real parts of the off diagonal elements of the scattering matrix. You will
notice that plot included demonstrates that the elements of the scattering
matrix are not smooth and a sign change occurs for certain energies (sign
changes for both the real and imaginary parts, plotting chosen for
simplicity).
Note that for the example I take a tiny range of energies, where the
conductance is the same and hence the number of active modes is constant as
a function of energy.
There may be a mistake in understanding on my part, however, I believe this
error is due to some sign error in the incoming or outgoing states in the
lead as this only effects individual modes (as you will see when you adjust
the system size, or move up to the next conductance plateau and the size of
the scattering matrix increases. My initial program was looking at graphene
for a large rectangular flake and I noticed this effect considerably). I
think this is a boundedness problem for certain functions which describes
the individual scattering states at certain energies.
Is this intentional behaviour for KWANT? As it means I can't perform
calculations such as calculating the Wigner-Smith time-delay matrix, among
other scattering characteristics.
Thank you for your time + best wishes,
Simon
My code:
#### KWANT tutorial code
from matplotlib import pyplot as plt
import kwant
import numpy as np
import scipy.linalg as la
# First, define the tight-binding system
sys = kwant.Builder()
# Here, we are only working with square lattices
a = 1
lat = kwant.lattice.square(a)
t = 1.0
W = 2
L = 2
# Define the scattering region
for i in xrange(L):
for j in xrange(W):
# On-site Hamiltonian
sys[lat(i, j)] = 4 * t
# Hopping in y-direction
if j > 0:
sys[lat(i, j), lat(i, j - 1)] = -t
# Hopping in x-direction
if i > 0:
sys[lat(i, j), lat(i - 1, j)] = -t
# Then, define and attach the leads:
# First the lead to the left
# (Note: TranslationalSymmetry takes a real-space vector)
sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
left_lead = kwant.Builder(sym_left_lead)
for j in xrange(W):
left_lead[lat(0, j)] = 4 * t
if j > 0:
left_lead[lat(0, j), lat(0, j - 1)] = -t
left_lead[lat(1, j), lat(0, j)] = -t
sys.attach_lead(left_lead)
# Then the lead to the right
sym_right_lead = kwant.TranslationalSymmetry((a, 0))
right_lead = kwant.Builder(sym_right_lead)
for j in xrange(W):
right_lead[lat(0, j)] = 4 * t
if j > 0:
right_lead[lat(0, j), lat(0, j - 1)] = -t
right_lead[lat(1, j), lat(0, j)] = -t
sys.attach_lead(right_lead)
# Plot it, to make sure it's OK
kwant.plot(sys)
# Finalize the system
sys = sys.finalized()
# Now that we have the system, we can compute conductance
emin, emax = 1.05, 1.06#3.00
points = 21
energies = np.linspace(emin,emax,points)
element_0_1, element_1_0 = [], []
for ie in energies:
energy = ie
# I only consider element 0,1 and 1,0 as 0,0 and 1,1 are vanishingly
small and look like floating point errors.
#This gets the smartix and looks at the real parts of the off-diagonal
elements
smatrix = kwant.smatrix(sys, energy).data
element_0_1.append(smatrix[0,1].real)
element_1_0.append(smatrix[1,0].real)
print smatrix[0,1].real
print smartix[1,0].real
fig2 = plt.figure()
ax2 = fig2.add_subplot(1,1,1)
ax2.plot(element_0_1, 'ro')
ax2.plot(element_1_0, 'bo')
plt.show()