Dear Kwant users,
I would like to plot the band structure of my system under a magnetic
field. Based on the code demostrated in
http://topocondmat.org/w3_pump_QHE/Laughlinargument.html, I plot the band
structure but I cannot observe any Landau level like spectrum. I suspect it
has something to do with the line
p = SimpleNamespace(t=1.0, mu=0.3, mu_lead=0.3, B=0.5)
bands = kwant.physics.Bands(flead,args = [p])
Can the B-field be passed to the Hamiltonian in this way? The complete code
is attached below. Thanks for your help.
Best,
Johnny
#------My code-----------#
from __future__ import division
import numpy as np
import kwant
from matplotlib import pyplot as plt
from math import pi
class SimpleNamespace(object):
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
def qhe_corbino(r_out=100, r_in=65, w_lead=10):
"""Create corbino disk.
Square lattice, one orbital per site.
Returns kwant system.
Arguments required in onsite/hoppings:
t, mu, mu_lead, B, phi
"""
# ring shape
def ring(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return (r_in ** 2 < rsq < r_out ** 2)
# Onsite and hoppings
def onsite(site, p):
(x, y) = site.pos
return 4 * p.t - p.mu
def crosses_branchcut(hop):
x1, y1 = hop[0].pos
x2, y2 = hop[1].pos
return y1 < 0 and x1 > 0.5 and x2 < 0.5
def hopping(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
# Check for correctness!
return -p.t * np.exp(-0.5j * p.B * (x1 - x2) * (y1 + y2))
def branchcut_hopping(site1, site2, p):
return hopping(site1, site2, p) * np.exp(1j * p.phi)
# Building system
lat = kwant.lattice.square()
sys = kwant.Builder()
sys[lat.shape(ring, (0, r_in + 1))] = onsite
sys[lat.neighbors()] = hopping
# adding special hoppings
def hops_across_cut(sys):
for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(sys):
if crosses_branchcut(hop):
yield hop
sys[hops_across_cut] = branchcut_hopping
# Attaching leads
sym_lead = kwant.TranslationalSymmetry((-1, 0))
lead = kwant.Builder(sym_lead)
def lead_shape(pos):
(x, y) = pos
return (-w_lead / 2 < y < w_lead / 2)
lead[lat.shape(lead_shape, (0, 0))] = lambda site, p: 4 * p.t -
p.mu_lead
lead[lat.neighbors()] = lambda site1, site2, p: -p.t
#### Attach the leads and return the system. ####
sys.attach_lead(lead)
sys.attach_lead(lead, origin=lat(0, 0))
return sys
def total_charge(value_array):
determinants = [np.linalg.det(s) for s in value_array]
charge = np.cumsum(np.angle(determinants / np.roll(determinants, 1)))
return charge / (2 * np.pi)
def qhe_ribbon(W, periodic=False):
""" Creates ribbon system
If we have periodic boundary conditions, the flux through a single
unit cell is quantized.
"""
W = 2 * (W // 2)
def ribbon_shape(pos):
(x, y) = pos
return (-W / 2 <= y <= W / 2)
def onsite(site, p):
(x, y) = site.pos
return 4 * p.t - p.mu
def hopping(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
return -p.t * np.exp(-0.5j * p.B * (x1 - x2) * (y1 + y2))
def hopping_periodic(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
return -p.t * np.exp(-1j * np.pi / (W + 1) * np.round((W + 1) * p.B
/ (2 * np.pi)) * (x1 - x2) * (y1 + y2))
lat = kwant.lattice.square()
sym_sys = kwant.TranslationalSymmetry((-1, 0))
sys = kwant.Builder(sym_sys)
sys[lat.shape(ribbon_shape, (0, 0))] = onsite
if periodic:
sys[lat.neighbors()] = hopping_periodic
sys[lat(0, - W / 2), lat(0, + W / 2)] = hopping_periodic
else:
sys[lat.neighbors()] = hopping
return sys
# Quantum hall bar codes
def qhe_hall_bar(L=50, W=10, w_lead=10, w_vert_lead=None):
"""Create a hall bar system.
Square lattice, one orbital per site.
Returns kwant system.
Arguments required in onsite/hoppings:
t, mu, mu_lead
"""
L = 2 * (L // 2)
W = 2 * (W // 2)
w_lead = 2 * (w_lead // 2)
if w_vert_lead is None:
w_vert_lead = w_lead
else:
w_vert_lead = 2 * (w_vert_lead // 2)
# bar shape
def bar(pos):
(x, y) = pos
return (x >= -L / 2 and x <= L / 2) and (y >= -W / 2 and y <= W / 2)
# Onsite and hoppings
def onsite(site, p):
(x, y) = site.pos
return 4 * p.t - p.mu
def hopping_Ax(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
return -p.t * np.exp(-0.5j * p.B * (x1 + x2) * (y1 - y2))
def make_lead_hop_y(x0):
def hopping_Ay(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
return -p.t * np.exp(-1j * p.B * x0 * (y1 - y2))
return hopping_Ay
def lead_hop_vert(site1, site2, p):
return -p.t
# Building system
lat = kwant.lattice.square()
sys = kwant.Builder()
sys[lat.shape(bar, (0, 0))] = onsite
sys[lat.neighbors()] = hopping_Ax
# Attaching leads
sym_lead = kwant.TranslationalSymmetry((-1, 0))
lead = kwant.Builder(sym_lead)
def lead_shape(pos):
(x, y) = pos
return (-w_lead / 2 <= y <= w_lead / 2)
lead_onsite = lambda site, p: 4 * p.t - p.mu_lead
sym_lead_vertical = kwant.TranslationalSymmetry((0, 1))
lead_vertical1 = kwant.Builder(sym_lead_vertical)
lead_vertical2 = kwant.Builder(sym_lead_vertical)
def lead_shape_vertical1(pos):
(x, y) = pos
return (-L / 4 - w_vert_lead / 2 <= x <= -L / 4 + w_vert_lead / 2)
def lead_shape_vertical2(pos):
(x, y) = pos
return (+L / 4 - w_vert_lead / 2 <= x <= +L / 4 + w_vert_lead / 2)
lead_vertical1[lat.shape(lead_shape_vertical1, (-L / 4, 0))] =
lead_onsite
lead_vertical1[lat.neighbors()] = lead_hop_vert
lead_vertical2[lat.shape(lead_shape_vertical2, (L / 4, 0))] =
lead_onsite
lead_vertical2[lat.neighbors()] = lead_hop_vert
sys.attach_lead(lead_vertical1)
sys.attach_lead(lead_vertical2)
sys.attach_lead(lead_vertical1.reversed())
sys.attach_lead(lead_vertical2.reversed())
lead[lat.shape(lead_shape, (-1, 0))] = lead_onsite
lead[lat.neighbors()] = make_lead_hop_y(-L / 2)
sys.attach_lead(lead)
lead = kwant.Builder(sym_lead)
lead[lat.shape(lead_shape, (-1, 0))] = lead_onsite
lead[lat.neighbors()] = make_lead_hop_y(L / 2)
sys.attach_lead(lead.reversed())
return sys
def calculate_sigmas(G):
# reduce by one dimension G -> G[temp, temp]
temp = [0, 1, 3, 4, 5]
G = G[temp, :]
G = G[:, temp]
# invert R = G^-1
# find out whether it is a numpy object
r = np.linalg.inv(G)
# Voltages follow: V = R I[temp]
V = r.dot(np.array([0, 0, 0, 1, -1]))
# Completely solved the six terminal system.
# Consider the 2x2 conductance now: Use I = sigma U
E_x = V[1] - V[0]
E_y = V[1] - V[3]
# formula above
sigma_xx = E_x / (E_x**2 + E_y**2)
sigma_xy = E_y / (E_x**2 + E_y**2)
return sigma_xx, sigma_xy
def plot_bandstructure(flead, momenta):
bands = kwant.physics.Bands(flead,args = [p])
energies = [bands(k) for k in momenta]
plt.figure()
plt.plot(momenta, energies)
plt.xlabel(r"$k[a^{-1}]$")
plt.ylabel(r"$E [t]$")
plt.show()
sys = qhe_hall_bar(L=60, W=100, w_lead=90, w_vert_lead=28).finalized()
kwant.plot(sys)
p = SimpleNamespace(t=1.0, mu=0.3, mu_lead=0.3, B=0.5)
Bs = np.linspace(0.02, 0.15, 200)
num_leads = len(sys.leads)
momenta = [-pi + 0.02 * pi * i for i in range(101)]
plot_bandstructure(sys.leads[2], momenta)
