Dear Pablo, thank You for Your answer - it clarifies a lot, but at the same time a few questions in my head appear. (At the beginning I'll mention that I'm more familiar with analytical calculation thus maybe my questions are trivial for the people familiar with numerical methodes, but I just want to understand some issues and figure out what is going in the results which I've obtained here.) 1. Now I have noticed that in Your previous reply You had changed the hoppings in the system in this way: <code>-t*sigma_0</code> -> <code>-t*sigma_z</code> and without this change the energy gap does not appear. I don't understand this change, because for me it creates a different system: this part originates from the kinetic term of the Hamiltonian and this term is defined with the identity matrix (sigma_0). Could You explain me why we can introduce that change? 2. Regarding the noises in xy-conductivity outside the gap in the finite system. You wrote about the tricks how to reduce this noise by changing the parametres of computation in KPM method - I checked, the noise is in fact smaller. But I want to find out how the xy-component will behave in the infinite system using the scattering-matrix method (will I get the results expected for the infinite system?). Thus, I've added the top and bottom leads and with the <code>smatrix.transmission()</code> method I've calulated the xx- and xy-conductance (based on 2.3.1 tutorial). Unfortunately, the results are unclear for me, because: a) adding the top and bottom leads to the system cancels the quantization in xx-conductance - I do not understand why adding leads perpendicular to the current flow does have the influence on the longitudinal conductance, b) the energy dependence of the both, the xx- and xy-conductance, increasing with oscillations - I don't understand this behaviour either, c) if I replace the hopping terms as You recommended, means <code>-t*sigma_0</code> -> <code>-t*sigma_z</code>, the results are unexpected either (0 in conductance till some energy, then random positive values). 3. The third thing, where I have a problem, is about the band structure. Taking a pen and a piece of paper: if I defined the system describing by the Hamiltonian which is the matrix 2x2, I will get two eigenvalues corresponding to two branches of the dispersion relation. Here, by calculating the band structure of the lead (describing by 2x2 matrix written out as on-site and hopping parameters), I get a lot of bands (which number increases with the size of the lead/system). I suppose that in this notation every unit cell/site produces a vector, and thus eigenvalue, and in that way we obtain so plentiful band structure, but it's not clear for me. Below I attach the code for the 2. point. Best regards, Anna ----------------------------------------------------- <code> import kwant from matplotlib import pyplot import tinyarray # we define the identity matrix and Pauli matrices 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]]) def make_system(t=1.0, alpha=1, e_z=0.5, W=30, L=30, a=1): # for simplicity we set 'a' equal to 1 lat = kwant.lattice.square() syst = kwant.Builder() # Zeeman energy adds to the on-site term syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z # hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) #syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+ 1j*alpha*sigma_y/(2*a) #change from the Pablo's code # hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) #syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z + 1j*alpha*sigma_x/(2*a) #change from the Pablo's code lead = kwant.Builder(kwant.TranslationalSymmetry((-a,0))) lead[(lat(0,j) for j in range(W))] = 4*t*sigma_0 + e_z*sigma_z lead[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction lead[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction # attach the left and right leads syst.attach_lead(lead) # left syst.attach_lead(lead.reversed()) # right # create the leads (top and bottom) lead2 = kwant.Builder(kwant.TranslationalSymmetry((0,a))) lead2[(lat(i,0) for i in range(L))] = 4*t*sigma_0 + e_z*sigma_z lead2[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction lead2[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction # attach the leads (top and bottom) syst.attach_lead(lead2) # top - ADDING THIS ELECTRODE changes the xx-component of conductivity tensor syst.attach_lead(lead2.reversed()) # bottom # return the finalized system return syst def plot_conductance(syst, energies, firstLead, secondLead): data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(secondLead,firstLead)) # form lead 'firstLead' to 'secondLead' pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel(r"conductance [e$^2$/h]") pyplot.title("form " + str(firstLead) + " lead to " + str(secondLead) + " lead") pyplot.show() def main(): syst = make_system() kwant.plot(syst); # check that the system looks as intended syst = syst.finalized() # calculate the band structure (of the lead) # the results are the same for each (left,right,top,bottom) lead kwant.plotter.bands(syst.leads[0]) print('leads:\t 0 - left \n\t 1 - right \n\t 2 - top \n\t 3 - bottom') # longitudinal conductivity (right-left = top-bottom) plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=0, secondLead=1) # from right to left #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=2, secondLead=3) # from bottom to top # transverse conductivity (reversible values, means top-left = left-top) plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=0, secondLead=2) # form top to left #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=2, secondLead=0) # form left to top #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=1, secondLead=2) # form top to right # standard Python construct to execute 'main' if the program is used as a script if __name__ == '__main__': main() </code>