Hello Kwant community, 
I want to calculate the local Chern marker in kwant, roughly a real-space 
approach to topological numbers such as the Chern number. Basically, one needs 
to calculate the expectation value of the commutator [x, y] of position 
operators projected to the occupied states. The main reference is Eqs.(6-9) in 
[ http://dx.doi.org/10.1103/PhysRevB.84.241106 ], which has been applied in 
many other cases.

My understanding is that kwant can define x,y as ???density operators??? and 
then this task should fit into the kwant.kpm.correlator framework, because it 
is a correlation function at zero temperature between the two operators. So I 
tried in the following code a few ways (some commented) for the simplest 
square-lattice quantum Hall effect (QHE). But it does not seem to give a sign 
of quantization, although I'm sure that the same system gives reasonable QHE in 
conductivity and conductance calculation. E.g., for 1/B in [3, 12], it should 
show the C=1, C=2 plateaus. I feel that I probably missed or misunderstood 
something basic here. Is my use of operator and kpm.correlator correct here? 
Thank you for your help!

from cmath import exp
import numpy as np
import matplotlib.pyplot as plt
import kwant
import scipy.sparse.linalg as sla
import scipy.sparse as sparse
import tinyarray

Nx, Ny = 101, 101
Lx, Ly = Nx-1, Ny-1
edge=18
eps = 1e-7
L_sys = [Lx, Ly] 

Ef=0.4

def make_system(a=1, t=1):
    lat = kwant.lattice.square(a, norbs=1)
    syst = kwant.Builder()
    def fluxx(site1, site2, Bvec):
        pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)]
        return exp(-1j * Bvec * pos[1] )
    def fluxy(site1, site2, Bvec):
        pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)]
        return 1 #exp(-1j * (-Bvec * pos[0]))

    def hopx(site1, site2, Bvec):
        # The magnetic field is controlled by the parameter Bvec
        return fluxx(site1, site2, Bvec) * t
    def hopy(site1, site2, Bvec):
        # The magnetic field is controlled by the parameter Bvec
        return fluxy(site1, site2, Bvec) * t
    def onsite(site):
        return 4*t

    #### Define the scattering region. ####
    syst[(lat(x, y) for x in range(L_sys[0]) for y in range(L_sys[1]))] = onsite
    # hoppings in x-direction
    syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
    # hoppings in y-directions
    syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy

    # Finalize the system.
    fsyst = syst.finalized()
    return lat, fsyst

lat, fsyst = make_system()

sites = fsyst.sites
X = sparse.diags([site.pos[0] for site in sites])
Y = sparse.diags([site.pos[1] for site in sites])


def where(site):
    pos=site.pos
    x0_min, y0_min = edge-eps, edge-eps
    x0_max, y0_max= Lx-edge+eps, Ly-edge+eps
    return x0_min < pos[0] < x0_max and y0_min < pos[1] < y0_max

def op0(site):
    return site.pos[0]

def op1(site):
    return site.pos[1]

def chern(lat, fsyst):
    # correlator=kwant.kpm.Correlator(fsyst, params=params, 
operator1=kwant.operator.Density(fsyst, onsite=op0, where=where),\
    #                      operator2=kwant.operator.Density(fsyst, onsite=op1, 
where=where),\
    #                      num_vectors=10, num_moments=400, 
vector_factory=kwant.kpm.RandomVectors(fsyst, where=where),\
    #                      bounds=(0, 8.), eps=0.01, rng=0, kernel=None, 
mean=True, accumulate_vectors=False)
    
    
    correlator=kwant.kpm.Correlator(fsyst, params=params, 
operator1=kwant.operator.Density(fsyst, onsite=op0),\
                         operator2=kwant.operator.Density(fsyst, onsite=op1),\
                         num_vectors=10, num_moments=400, 
vector_factory=kwant.kpm.RandomVectors(fsyst, where=where),\
                         bounds=(0, 8.), eps=0.01, rng=0, mean=True, 
accumulate_vectors=False)
    
    
    # correlator=kwant.kpm.Correlator(fsyst, params=params, operator1=X,\
    #                      operator2=Y,\
    #                      num_vectors=10, num_moments=400, 
vector_factory=kwant.kpm.RandomVectors(fsyst, where=where),\
    #                      bounds=(0, 8.), eps=0.01, rng=0, kernel=None, 
mean=True, accumulate_vectors=False)

    c = correlator(mu=Ef, temperature=0.0) / ((Lx-2*edge)*(Ly-2*edge))
    print(2*np.pi*c)
    print("Chern number: ", np.real(np.pi*1j*(c-c.conj()))) # Chern number
    
Bfields = []
for BInv in np.arange(3., 12., 2):
    Bfields.append(1/BInv)
NB = len(Bfields)
params = dict(Bvec=0)
for B in Bfields:
    params['Bvec'] = B
    print("B: ", B)
    chern(lat, fsyst)

Reply via email to