Hello all, I'm trying to implement the BHZ model already in the tight-biding model according to the paper [Phys. Rev. B 96, 155302 (2017) or arXiv:1707.03624]. I use the equations (5),(6) and (7), however with no magnetic field. I'm using Kwant's documentation example 2.10 as a guideline (For what I have to do next it's more convenient that I don't use the discretize method). After some struggling I got the band structure of both mine and Kwant's example to be the same, however my spin density is reversed, it's up where it's supposed to be down. So far I haven't figured out why. Also, the main difference between the implementation is only the make_syst part, all the plotting in my code are snippets from Kwant's discretize.py; My code is below: import kwant import matplotlib.pyplot as plt import numpy as np import random as rd from math import * import matplotlib as mpl import tinyarray 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]]) Id = np.eye(4) G1 = np.kron(sigma_z,sigma_x) G2 = -np.kron(sigma_0,sigma_y) G5 = np.kron(sigma_0,sigma_z) def make_syst(a=20,A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0,L=200,W=200): lat = kwant.lattice.square(a,norbs=4) syst = kwant.Builder() def square(pos): (x, y) = pos return (0 <= y < W and 0 <= x < L) def onsite(site): (x, y) = site.pos return C*Id+M*G5-(4/a**2)*(D*Id+B*G5) #equation 5 def hoppingx(site1,site2): return (1/a**2)*(D*Id+B*G5)+(1j/(2*a))*(A*G1) #eq 6 def hoppingy(site1,site2): return (1/a**2)*(D*Id+B*G5)-(1j/(2*a))*(A*G2) #eq 7 syst[lat.shape(square,(0,0))] = onsite syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hoppingx syst[kwant.builder.HoppingKind((0,1), lat, lat)] = hoppingy def lead_shape(pos): (x, y) = pos return (0 <= y < W) #### Define the left lead. #### lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0))) lead[lat.shape(lead_shape,(0,0))] = onsite # hoppings in x-direction lead[kwant.builder.HoppingKind((1, 0), lat, lat)] =hoppingx # hoppings in y-directions lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = hoppingy return syst,[lead,lead.reversed()] def main(): syst,leads = make_syst(W=1000,L=2000) # Attach the leads to the system. for lead in leads: syst.attach_lead(lead) # Finalize the system. syst = syst.finalized() # Check that the system looks as intended. #kwant.plot(syst) # Compute the band structure of lead 0. kwant.plotter.bands(syst.leads[0], momenta=np.linspace(-0.3, 0.3, 201), show=False) plt.grid() plt.xlim(-.3, 0.3) plt.ylim(-0.05, 0.05) plt.xlabel('momentum [1/A]') plt.ylabel('energy [eV]') plt.show() wf = kwant.wave_function(syst, energy=0) # prepare density operators sigma_z = np.array([[1, 0], [0, -1]]) prob_density = kwant.operator.Density(syst, np.eye(4)) spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2))) # calculate expectation values and plot them wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(16, 4)) kwant.plotter.map(syst, wf_sqr, ax=ax1) kwant.plotter.map(syst, rho_sz, ax=ax2) ax = ax1 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax = ax2 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax1.set_title('Probability density') ax2.set_title('Spin density') plt.show() # Call the main function if the script gets executed (as opposed to imported). # See http://docs.python.org/library/__main__.html. if __name__ == '__main__': main() Regards, Victor Regis