Hi All,
I am trying to reproduce Fig. 4 in PRL 97, 067007 (2006) - Specular Andreev
Reflection in Graphene. This had been done in KWANT in the PhD thesis of Tibor
Sekera - "Quantum transport of fermions in honeycomb lattices and cold atomic
systems", chapter 4. My code is the following. Can Anyone help me? I am not
getting same result. I am trying the same code this question has been asked
before here
(https://mail.python.org/archives/list/kwant-discuss@python.org/thread/3JH5W2YUYPVPDVTAE6I57JFLG4UDAZN5/),
but no fruitful outcome was gained. I shall be highly grateful if anyone can
help me.
################################################################
import kwant, tinyarray
from math import pi, sqrt
from matplotlib import pyplot as plt
import numpy as np
plt.rcParams['figure.dpi']=150
####___Fundamental 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]])
#Fix hopping, Fermi energy, and superconducting gap
t = 1
# Ef = 0
Ef_S=2
Delta= 2
a=1
V0=0.5
#Length and Width of scattering region (the superconductor)
L = 20
W = 10
def make_system():
lat =
kwant.lattice.general([(sqrt(3)*0.5,0.5),(0,1)],[(0,0),(1/(2*sqrt(3)),1/2)],norbs
= 4)
a, b = 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_S * -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): # left lead
x, y = pos
return -W <=y <=W
sym_up = kwant.TranslationalSymmetry(lat.vec((-2,1))) # Symmetry along the
x-axis
# 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_S * -1 * np.kron(tau_z, s_0)
lead0[lat.neighbors()] = t * -1 * np.kron(tau_z, s_0)
########################################################################################################################
def lead1_shape(pos): # right
x, y = pos
return -W <=y <=W
sym_down = kwant.TranslationalSymmetry(lat.vec((2, -1)))
# 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 * np.kron(tau_z, s_0)
lead1[lat.neighbors()] = t * -1 * np.kron(tau_z, s_0)
return syst, [lead0, lead1]
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)))
plt.figure()
plt.plot(energies, data,'*-r')
plt.xlabel("energy [t]")
plt.ylabel("conductance [e^2/h]")
plt.show()
#Plot the band in the lead
#Make my system and attach the leads
def main():
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.
kwant.plot(system,fig_size=(5,5))
system = system.finalized()
# #Evaluate the band in the lead
L_1 = leads[0].finalized()
# bands = kwant.physics.Bands(L_1)
# #Plot the band in the lead
# ks = np.linspace(-pi,pi,num=101)
# energies=[bands(k) for k in ks]
# f, ax = plt.subplots()
# ax.plot(ks,energies)
# plt.show()
plot_conductance(system, energies=[0.01 * i+0.00001 for i in range(200)])
main()
