Hi all,
I am somewhat new to Kwant and am trying to achieve the following:
I have a 2D 3-band continuum model that I want to see edge states of using band
and density diagrams via a 1D real-space strip. I make the system finite in the
x direction (implying translational symmetry in the y direction, making k_y a
good quantum number and so constant). I setup the domain wall by choosing
between two values for a potential, dependent on whether the site is before
length/2 or after. I use kwant.continuum.discretize() for builder().
However, I am having trouble with the setup in general. For one, I cannot seem
to apply translational symmetry to builder() because I am not sure what vector
to use (I tried (0,1) and (1,0)), but to no avail. Second, if I proceed without
translational symmetry, kwant.plotter.bands(syst, show=False) gives me
"TypeError: Expecting an instance of InfiniteSystem." When I browsed the email
archives, people suggested syst.leads[0], etc: but the issue is that I don't
have these leads explicitly and that gives me an error. My understanding is
that I can get a proper band diagram showing edge states if I plot eigenvalues
from looping through different values of k_y that I fixed earlier. However, I
don't think I saw examples of the Kwant examples doing it this way.
So, it seems as if I have confused some basic applications in Kwant. Would
someone mind clarifying my confusions for me? All I want to see is the a
typical band structure diagram for this 3-band continuum model in a domain wall
setup. Following is my code (I know that the given Hamiltonian will not give
edge states, but I just want to get the code working for starters). Thank you
for your time!
# Minimal example
import kwant.continuum
import kwant
import numpy as np
import matplotlib.pyplot as plt
import tinyarray
import scipy
lat_const = 3.19 # lattice constant
hami = """
[[2* t0 *
((1-(2*((1/2)*k_x*a))**2/2)+2*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2))+e1,
M1(x,y)-2*sqrt(3)*t2*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2)) +
2*1j*t1*((2*((1/2)*k_x*a)) +
(((1/2)*k_x*a))*(1-((momentumY*a*sqrt(3)/2))**2/2)),
2*t2*((1-(2*((1/2)*k_x*a))**2/2) -
(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)) +
2*sqrt(3)*1j*t1*(1-(((1/2)*k_x*a))**2/2)*((momentumY*a*sqrt(3)/2))],
[-M1(x,y)-2*sqrt(3)*t2*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2)) -
2*1j*t1*((2*((1/2)*k_x*a)) +
(((1/2)*k_x*a))*(1-((momentumY*a*sqrt(3)/2))**2/2)),
2*t11*(1-(2*((1/2)*k_x*a))**2/2) +
(t11+3*t22)*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)+e2,
-M2(x,y)+sqrt(3)*(t22-t11)*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2))+4*1j*t12*(((1/2)*k_x*a))*((1-(((1/2)*k_x*a))**2/2)-(1-((momentumY*a*sqrt(3)/2))**2/2))],
[2*t2*((1-(2*((1/2)*k_x*a))**2/2) -
(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)) -
2*sqrt(3)*1j*t1*(1-(((1/2)*k_x*a))**2/2)*((momentumY*a*sqrt(3)/2)),
M2(x,y)+sqrt(3)*(t22-t11)*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2))-4*1j*t12*(((1/2)*k_x*a))*((1-(((1/2)*k_x*a))**2/2)-(1-((momentumY*a*sqrt(3)/2))**2/2)),
2*t22*(1-(2*((1/2)*k_x*a))**2/2) + (3*t11+
t22)*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)+e2]]
"""
template = kwant.continuum.discretize(hami, grid = lat_const)
#print(template)
L = 100
def wire_shape(site):
(x, ) = site.pos
return (0<=x<L)
val1 = 10.07822923*1j
val2 = 12.42088476*1j
def pot1(x,y=0):
return[0,val1][x<L/2]
def pot2(x,y=0):
return[0,val2][x<L/2]
syst = kwant.Builder();
syst.fill(template, wire_shape,(0,));
syst = syst.finalized();
def get_ham(momY):
ham = syst.hamiltonian_submatrix(params=dict(M1=pot1,M2=pot2, momentumY =
momY, y=0, a = 3.19, e1 = 1.046, e2 = 2.104, t0 = -0.184, t1 = 0.401, t2 =
0.507, t11 = 0.218, t12 = 0.338, t22 = 0.057),sparse = False)
return ham
#plt.spy(ham) #block off-diagonal
#plt.show()
# Band structure diagram???
# different momentumY's to iterate over
v1=-2*np.pi; v2=2*np.pi; v_segs=10;
vary_range = np.linspace(v1, v2, num=v_segs) # np.linspace(kx0_1, kx0_2,
num=kx_segs, retstep=True)
plt.figure()
for i in range(v_segs):
ham = get_ham(vary_range[i])
ev,evecs = scipy.linalg.eigh(ham)
plt.plot(ev) # plot eigenvalues
plt.show()
# Wavefunction density diagram??? Should show peak at edge states, if there are
any.
plt.figure()
for i in range(v_segs):
ham = get_ham(vary_range[i])
ev,evecs = scipy.linalg.eigh(ham)
wavefunction = evecs[:, len(ham)-1]
up,middle,down = wavefunction[::3],wavefunction[1::3],wavefunction[2::3]
density = np.abs(up)**2 + np.abs(middle)**2 + np.abs(down)**2
plt.plot(density)
##kwant.plotter.map(syst, density, show=False)
plt.show()
Best,
Tharindu
