Hello Kwant community,
I have played with a Weyl semimetal system for a while, using Kwant's conductance features (with leads). And I do get the quantized Hall response under magnetic field as reported in this work [PRL 125, 036602].
Now with exactly the same system (size, Fermi energy, hopping, etc.), I remove the leads and try the kpm.conductivity features. However, there is no sign of quantization at all, as far as I've tried with LocalVectors or RandomVectors and tuning num_moments or so. To make sure of my use of kpm.conductivity, I tried another known unquantized Hall effect (sigma_xy) in this system without magnetic field and it looked not unreasonable; I also changed the system to a simple 2D QHE one and the same code worked well to see quantization.
This seems rather unclear to me. Could you please help with how should I correctly calculate the conductivity in Kwant here? The code below plots sigma_xz and sigma_zx with respect to a few inverse magnetic fields.
Thanks!
from cmath import exp
from math import cos, sin, pi
import numpy as np
import time
import sys
import kwant
# For plotting
from matplotlib import pyplot as plt
# For matrix support
import tinyarray
# define Pauli-matrices for convenience
sigma_0 = tinyarray.array([[1, 0], [0, 1]])
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])
Lx=25
Ly=20
Lz=25
D=0.02
D1=D
D2=4*D
A=1
M=0.3
kw=pi/2.0
E_F=2*D2*(1-cos(kw)) # Fermi energy
norbs=2
def make_system(a=1, t=1):
lat = kwant.lattice.cubic(a, norbs=norbs)
syst = kwant.Builder()
def fluxx(site1, site2, B):
pos = [(site1.pos[i]+site2.pos[i])/2.0 for i in range(3)]
return exp(-1j * (B * pos[2]))
def hopx(site1, site2, B):
return fluxx(site1, site2, B) * (-D2 * sigma_0 + 1j * A/2 * sigma_x + M * sigma_z)
def hopy(site1, site2):
return (-D1 * sigma_0 + 1j * A/2 * sigma_y + M * sigma_z)
def hopz(site1, site2):
return (-D2 * sigma_0 + M * sigma_z)
def onsite(site):
return (2 * D1 + 2 * D2 * 2) * sigma_0 + 2 * M * ((1 - cos(kw)) - 3) * sigma_z
syst[(lat(x, y, z) for x in range(Lx) for y in range(Ly) for z in range(Lz))] = onsite
# hoppings in x-direction
syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hopx
# hoppings in y-directions
syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hopy
# hoppings in z-directions
syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hopz
return lat, syst.finalized()
def plot_curves(labels_to_data):
plt.figure(figsize=(12,10))
for label, (x, y) in labels_to_data:
plt.plot(x, y, label=label, linewidth=2)
plt.legend(loc=2, framealpha=0.5)
plt.xlabel("1/B")
plt.ylabel(r'$\sigma$')
plt.show()
plt.clf()
def cond(lat, syst, random_vecs_flag=False):
center_width = min(Lx,Ly,Lz)//20
center_pos = [Lx//2,Ly//2,Lz//2]
def center(s):
return (np.abs(s.pos[0]-center_pos[0]) < center_width) and (np.abs(s.pos[1]-center_pos[1]) < center_width) and (np.abs(s.pos[2]-center_pos[2]) < center_width)
num_mmts = 200;
if random_vecs_flag: #use RandomVecs
center_vecs = kwant.kpm.RandomVectors(syst, where=center)
one_vec = next(center_vecs)
num_vecs = 20
else: # use LocalVecs
center_vecs = np.array(list(kwant.kpm.LocalVectors(syst, where=center)))
one_vec = center_vecs[0]
num_vecs = len(center_vecs)
area_per_orb = np.abs(np.linalg.det(lat.prim_vecs))
norm = area_per_orb * np.linalg.norm(one_vec) ** 2 / norbs
Bfields = []
Binvs = np.arange(2.4, 65.2, 30.5)
for BInv in Binvs:
Bfields.append(1/BInv)
params = dict(B=0)
cond_array_zx = []
cond_array_xz = []
for B in Bfields:
params['B'] = B
cond_zx = kwant.kpm.conductivity(syst, params=params, alpha='z', beta='x', mean=True, rng=0, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs)
cond_xz = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='z', mean=True, rng=0, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs)
cond_array_zx.append(cond_zx(mu=E_F, temperature=0.0)/ norm)
cond_array_xz.append(cond_xz(mu=E_F, temperature=0.0)/ norm)
cond_array_zx = np.array(cond_array_zx) * Ly # 3D conductivity * Ly gives the sheet 2D-like conductivity supposed to be quantized
cond_array_xz = np.array(cond_array_xz) * Ly
print(cond_array_zx)
print(cond_array_xz)
print('Calculating duration (s) until finalising system: ', time.time() - global_start_time)
plot_curves([(r'$\sigma_{zx}$', (Binvs, cond_array_zx)), (r'$\sigma_{xz}$', (Binvs, cond_array_xz))])
def main():
lat, syst = make_system()
# Check that the system looks as intended.
#kwant.plot(syst, file="lattice_fig.pdf")
cond(lat, syst)
if __name__ == '__main__':
global_start_time = time.time()
main()