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