Hello everyone, I am trying to implement Kitaev model in Kwant via the Hamiltonian: H = -\mu \sum_{n=1}^{N} c_{n}^{\dagger}c_{n} - t \sum_{n=1}^{N-1}(c_{n+1}^{\dagger}c_n + h.c.) + \Delta \sum_{n=1}^{N-1}(c_{n}c_{n+1} + h.c.) The first two terms are simple, but I am stucked in the superconducting p-wave pairing term. I've tried to implement it by doing sys[((lat_e(x+1, 0), lat_h(x, 0)) for x in range(0,L-1))] = Delta which seems somehow correct once I got the zero modes at the edges for the Kitaev limit \mu=0 and t=\Delta. But once I vary \mu I lost one of the edge states, and also it does not close the gap, which makes me belive that I am dealing with a trivial phase. Does anyone ever done it? Can you see something wrong with my approach? Cheers, Antônio. The complete code is the following: ==================================================== from cmath import exp import numpy as np import kwant # For eigenvalue computation import scipy.sparse.linalg as sla # For plotting from matplotlib import pyplot def make_system(a=1.0, W=1, L=50, mu=1.5, Delta=1.0, t=1.0): # Start with an empty tight-binding system and two square lattices, # corresponding to electron and hole degree of freedom lat_e = kwant.lattice.square(a=1.0, name='e') lat_h = kwant.lattice.square(a=1.0, name='h') sys = kwant.Builder() #### Define the scattering region. #### sys[(lat_e(x, y) for x in range(L) for y in range(W))] = - mu sys[(lat_h(x, y) for x in range(L) for y in range(W))] = mu # hoppings for both electrons and holes sys[lat_e.neighbors()] = -t sys[lat_h.neighbors()] = t # Superconducting order parameter enters as hopping between # electrons and holes sys[((lat_e(x+1, y), lat_h(x, y)) for x in range(0,L-1) for y in range(W))] = Delta #### Define the leads. #### # Symmetry for the left leads. sym_left = kwant.TranslationalSymmetry((-a, 0)) # left electron lead lead0 = kwant.Builder(sym_left) lead0[(lat_e(0, j) for j in range(W))] = - mu lead0[lat_e.neighbors()] = t lead0[lat_h.neighbors()] = -t # left hole lead lead1 = kwant.Builder(sym_left) lead1[(lat_h(0, j) for j in range(W))] = mu lead1[lat_e.neighbors()] = t lead1[lat_h.neighbors()] = -t # Then the lead to the right # this one is superconducting and thus is comprised of electrons # AND holes sym_right = kwant.TranslationalSymmetry((a, 0)) lead2 = kwant.Builder(sym_right) lead2 += lead0 lead2 += lead1 lead2[((lat_e(1, j), lat_h(0, j)) for j in range(W))] = Delta #### Attach the leads and return the system. #### sys.attach_lead(lead2) sys.attach_lead(lead2.reversed()) return sys def plot_wave_function(sys, L=100): # Calculate the wave functions in the system. ham_mat = sys.hamiltonian_submatrix(sparse=True) evecs = sla.eigsh(ham_mat, k=50, which='SM')[1] # Plot the probability density of the 1st eigenmode. kwant.plotter.map(sys, np.abs(evecs[:, 0])**2, colorbar=False, a=1.0, oversampling=2) def main(): sys = make_system() # Check that the system looks as intended. kwant.plot(sys) # Finalize the system. sys = sys.finalized() plot_wave_function(sys) if __name__ == '__main__': main() =================================================================