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)