Hello All,
I have seen that this problem has been discussed in threads before and thanks
to them, I've managed to set up most of it. But I need help on one main issue
I'm seeing right now.
What is unusual right now is by setting any non-zero fermi-energy in the normal
metal lead, the bands split and a gap opens. I wonder if I am using the wrong
matrices.
Code below:
import kwant
import numpy as np
from matplotlib import pyplot as plt
from math import pi, sqrt, tanh
import tinyarray
% matplotlib notebook
####___Fundemental Constants___####
e = 1.60217662 *(10**-19) #
hbar = 1.0545718*(10**-34) #
h = 6.62607004*(10**-34) #
e_hb = e/hbar #
###################################
#Specifying 2x2 Pauli matrices acting on spin
s_x = tinyarray.array([[0,1],[1,0]])
s_y = tinyarray.array([[0,-1j],[1j,0]])
s_z = tinyarray.array([[1,0],[0,-1]])
s_0 = tinyarray.array([[1,0],[0,1]])
#Specifying 2x2 Pauli matrices acting on particle-hole
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]])
delta_gap = 0.025
class SimpleNamespace(object):
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
# Ef_S is the fermi level in the superconducting region
# Ef is the fermi level in the normal graphene region
def make_system(a=0.142, W=1, L=1, t=2.71, tp=-0.1, V0=0.001, Ef=delta_gap*1,
Ef_S=delta_gap*10, Delta=delta_gap):
delta_t = lambda theta: Delta * np.cos(theta)
lat = kwant.lattice.general([(a * 1, 0), (a * 0.5, a * 0.5 * np.sqrt(3))],
[(0, 0), (0, a * 1 / np.sqrt(3))],
norbs=4)
a1, b1 = lat.sublattices
def channel(pos):
x, y = pos
return -W <= y <= W and -L <= x <= L
syst = kwant.Builder()
syst[lat.shape(channel, (0, 0))] = Ef * -1 * np.kron(tau_z,s_0)
syst[lat.neighbors()] = t *-1* np.kron(tau_z,s_0)
syst.eradicate_dangling()
def lead0_shape(pos): # top lead
x,y= pos
return -L <= x <= L
sym_up = kwant.TranslationalSymmetry(lat.vec((-1, 2)))
# Normal graphene lead-0
lead0 = kwant.Builder(sym_up, conservation_law= np.kron(tau_z,s_0),
particle_hole=np.kron(tau_x,s_0))
lead0[lat.shape(lead0_shape, (0, 0))] = Ef *-1* np.kron(tau_z,s_0)
lead0[lat.neighbors()] = t *-1* np.kron(tau_z,s_0)
def lead1_shape(pos): # bottom lead
x,y = pos
return -L <= x <= L
sym_down = kwant.TranslationalSymmetry(lat.vec((1, -2))) # bottom lead
# Superconducting graphene lead-1
lead1 = kwant.Builder(sym_down)
lead1[lat.shape(lead1_shape, (0, 0))] = (Ef_S + V0) *-1* np.kron(tau_z,s_0)
+ Delta * -1 * np.kron(tau_y,s_y)
lead1[lat.neighbors()] = t *-1* np.kron(tau_z,s_0)
return syst, [lead0,lead1]
def plot_bandstructure(flead, momenta, header):
bands = kwant.physics.Bands(flead)
energies = [bands(k) for k in momenta]
plt.figure(figsize=(9,6))
plt.plot(momenta, energies)
plt.xlabel("k $(1/a_0)$",fontsize=15)
plt.xlim(left=-1,right=1)
plt.ylabel("E (t)",fontsize=15)
plt.ylim(top=1,bottom=-1)
plt.tick_params(labelsize=12)
plt.title(header)
plt.show()
def plot_conductance(syst, energies):
# Compute transmission as a function of energy
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)))
plt.figure(figsize=(8,6))
plt.plot(energies/delta_gap, data)
plt.xlabel("$E/\Delta$",fontsize=15)
plt.ylabel("G ($e^2/h$)",fontsize=15)
# plt.ylim(bottom=0,top=2.5)
plt.tick_params(labelsize=15)
plt.show()
def main():
a=0.142
E = 0.001
system, leads = make_system()
# plot_system(system,a)
# Attach the leads to the system.
for lead in leads:
system.attach_lead(lead)
# Plot system with leads
def family_colors(site):
return 0 if site.family == a else 1
# kwant.plot(system, site_color=family_colors, site_lw=0.1, lead_site_lw=0,
colorbar=False)
# Finalize the system.
system = system.finalized()
# Compute the band structure of lead 0.
momenta = [-pi + 0.02 * pi * i for i in range(101)]
plot_bandstructure(system.leads[0], momenta,'Normal Lead')
plot_bandstructure(system.leads[1], momenta,'Superconducting Lead')
# Plot conductance.
energies = np.linspace(-0,2*delta_gap,201)
plot_conductance(system, energies)
if __name__ == '__main__':
main()
Any guidance or corrections will be appreciated!
Thank you,
Sharadh
