Hi Anton, Thanks for your help. I have resolved that issue myself and am writing about a current issue I am facing. *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? 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() On Wed, Sep 7, 2022 at 5:13 AM Anton Akhmerov <anton.akhmerov+kd@gmail.com> wrote:

Hi Binayyak,

Your description is insufficient to identify the source of your problem. Please see https://stackoverflow.com/help/how-to-ask and https://betatim.github.io/posts/getting-answers/ for guidance on how to ask technical questions.

Best, Anton

On Wed, 7 Sept 2022 at 11:11, Binayyak Roy <binayyr@g.clemson.edu> wrote:

...

Hey,

I have been trying to set up a 1D quantum wire and then calculate the

conductance. However when trying to calculate the smatrix it gives me the following error:

...

ValueError: Input needs to be a square matrix.

And my system that I am working with is defined as follows: lat =

kwant.lattice.chain(a)

...

How can I possibly resolve this problem? Or where should I be looking to

get an idea?

...

Sincerely, Binayyak