Hi, *I am trying to build a 2D quantum wire with both Rashba spin coupling and superconductivity*. I have been building on the code of superconductivity already mentioned in the kwant documentation pdf and added the rashba coupling terms to the code. However when I add the asymmetric terms of the rashba spin coupling it creates some sort of an error giving me an empty scattering matrix. What could be causing this? After going through the documentation I realize that it's probably got something to do with the fact that I have the wrong conservation law. I am not sure how to set that up? Sincerely, Binayyak I have attached my code below. import kwant import tinyarray import numpy as np # For plotting from matplotlib import pyplot tau_0 = tinyarray.array([[1, 0], [0, 1]]) tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]]) def make_system( a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4), mu=0.4, Delta=0.1, Deltapos=4, t=1.0, E_z=1e-1, # zeeman splitting field alpha=1e-3, # rashba coupling coefficient phs=True, ): lat = kwant.lattice.square(norbs=2) syst = kwant.Builder() #### Define the scattering region. #### # The superconducting order parameter couples electron and hole orbitals # on each site, and hence enters as an onsite potential. # The pairing is only included beyond the point 'Deltapos' in the scattering region. syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = ( 4 * t - mu + E_z ) * tau_z syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = ( 4 * t - mu + E_z ) * tau_z + Delta * tau_x # The tunnel barrier syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1]) for y in range(W))] = ( 4 * t + barrier - mu + E_z ) * tau_z # Hoppings syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t * tau_z + 1j * alpha * tau_y syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t * tau_z - 1j * alpha * tau_x #### Define the leads. #### # Left lead - normal, so the order parameter is zero. sym_left = kwant.TranslationalSymmetry((-a, 0)) # Specify the conservation law used to treat electrons and holes separately. # We only do this in the left lead, where the pairing is zero. lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y) lead0[(lat(0, j) for j in range(W))] = (4 * t - mu + E_z) * tau_z lead0[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t * tau_z + 1j * alpha * tau_y lead0[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t * tau_z - 1j * alpha * tau_x # Right lead - superconducting, so the order parameter is included. sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) lead1[(lat(0, j) for j in range(W))] = (4 * t - mu + E_z) * tau_z + Delta * tau_x lead1[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t * tau_z + 1j * alpha * tau_y lead1[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t * tau_z - 1j * alpha * tau_x #### Attach the leads and return the system. #### syst.attach_lead(lead0) syst.attach_lead(lead1) for i in syst.leads[0].interface: for j in syst.neighbors(i): syst[i, j] = 0.1 * tau_0 for i in syst.leads[1].interface: for j in syst.neighbors(i): syst[i, j] = 0.1 * tau_0 return syst def plot_conductance(syst, energies): # Compute conductance data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append( smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0)) ) pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel("conductance [e^2/h]") pyplot.show() def check_PHS(syst): # Scattering matrix s = kwant.smatrix(syst, energy=0) # Electron to electron block s_ee = s.submatrix((0, 0), (0, 0)) # Hole to hole block s_hh = s.submatrix((0, 1), (0, 1)) print("s_ee: \n", np.round(s_ee, 3)) print("s_hh: \n", np.round(s_hh[::-1, ::-1], 3)) print("s_ee - s_hh^*: \n", np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), "\n" ) # Electron to hole block s_he = s.submatrix((0, 1), (0, 0)) # Hole to electron block s_eh = s.submatrix((0, 0), (0, 1)) print("s_he: \n", np.round(s_he, 3)) print("s_eh: \n", np.round(s_eh[::-1, ::-1], 3)) print("s_he + s_eh^*: \n", np.round(s_he + s_eh[::-1, ::-1].conj(), 3)) def main(): syst = make_system(W=10) # Check that the system looks as intended. kwant.plot(syst) # Finalize the system. syst = syst.finalized() # Check particle-hole symmetry of the scattering matrix check_PHS(syst) # Compute and plot the conductance plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)]) # 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()