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()