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

Reply via email to