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()