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