Hi Simon,

Kwant algorithms do not adapt any particular convention in choosing the phases of the modes in the leads. That's why the scattering matrix at different energies can be discontinuous. You can access the wave functions of the modes through the resulting smatrix [1]. You can also easily modify the way Kwant calculates modes like we've done in [2].

Hope that helps, Anton

[1]: http://kwant-project.org/doc/1.0/reference/generated/kwant.solvers.common.SM... [2]: http://arxiv.org/abs/1408.1563 specifically http://arxiv.org/src/1408.1563v2/anc/trijunction.py

On Mon, Sep 28, 2015 at 7:41 PM, Simon Malzard smalzard01@gmail.com wrote:

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