Spin Hall effect in four terminal devices
From phys824
Jump to navigationJump to search
KWANT script below describes four-terminal graphene device, with gold adatoms in the central square, which generates spin Hall current in the transverse leads as a response to injected longitudinal charge current . The script output is the spin Hall angle defined as .
It may be required to install sympy package. For Windows users, download the .whl file from this link. Then, go to the directory with the .whl file, e.g.:
cd c:\Users\YOUR_USERNAME\Downloads
To install all the .whl-files in the current directory, execute
python -c "import pip, glob; pip.main(['install', '--no-deps'] + glob.glob('*.whl'))"
For more information, refer to Kwant installation instructions here.
from math import sqrt
import random
import itertools as it
import tinyarray as ta
import numpy as np
import matplotlib.pyplot as plt
import kwant
from sympy import S, Eq, solve, Symbol
class Honeycomb(kwant.lattice.Polyatomic):
"""Honeycomb lattice with methods for dealing with hexagons"""
def __init__(self, name=''):
prim_vecs = [[0.5, sqrt(3)/2], [1, 0]] # bravais lattice vectors
# offset the lattice so that it is symmetric around x and y axes
basis_vecs = [[-0.5, -1/sqrt(12)], [-0.5, 1/sqrt(12)]]
super(Honeycomb, self).__init__(prim_vecs, basis_vecs, name)
self.a, self.b = self.sublattices
def hexagon(self, tag):
""" Get sites belonging to hexagon with the given tag.
Returns sites in counter-clockwise order starting
from the lower-left site.
"""
tag = ta.array(tag)
# a-sites b-sites
deltas = [(0, 0), (0, 0),
(1, 0), (0, 1),
(0, 1), (-1, 1)]
lats = it.cycle(self.sublattices)
return (lat(*(tag + delta)) for lat, delta in zip(lats, deltas))
def hexagon_neighbors(self, tag, n=1):
""" Get n'th nearest neighbor hoppings within the hexagon with
the given tag.
"""
hex_sites = list(self.hexagon(tag))
return ((hex_sites[(i+n)%6], hex_sites[i%6]) for i in range(6))
def random_placement(builder, lattice, density):
""" Randomly selects hexagon tags where adatoms can be placed.
This avoids the edge case where adatoms would otherwise
be placed on incomplete hexagons at the system boundaries.
"""
assert 0 <= density <= 1
system_sites = builder.sites()
def hexagon_in_system(tag):
return all(site in system_sites for site in lattice.hexagon(tag))
# get set of tags for sites in system (i.e. tags from only one sublattice)
system_tags = (s.tag for s in system_sites if s.family == lattice.a)
# only allow tags which have complete hexagons in the system
valid_tags = list(filter(hexagon_in_system, system_tags))
N = int(density * len(valid_tags))
total_hexagons=len(valid_tags)
valuef=random.sample(valid_tags, N)
return valuef
def cross(W, L):
def shape(pos):
return ((-W <= pos[0] <= W and -L <= pos[1] <= L) or # vertical
(-W <= pos[1] <= W and -L <= pos[0] <= L)) # horizontal
return shape
## Pauli matrices ##
sigma_0 = ta.array([[1, 0], [0, 1]]) # identity
sigma_x = ta.array([[0, 1], [1, 0]])
sigma_y = ta.array([[0, -1j], [1j, 0]])
sigma_z = ta.array([[1, 0], [0, -1]])
## Hamiltonian value functions ##
def onsite_potential(site, params):
return params['ep'] * sigma_0
def potential_shift(site, params):
return params['mu'] * sigma_0
def kinetic(site_i, site_j, params):
return -params['gamma'] * sigma_0
def rashba(site_i, site_j, params):
d_ij = site_j.pos - site_i.pos
rashba = 1j * params['V_R'] * (sigma_x * d_ij[1] - sigma_y * d_ij[0])
return rashba + kinetic(site_i, site_j, params)
def spin_orbit(site_i, site_j, params):
so = (1j/3.) * params['V_I'] * sigma_z
return so
lat = Honeycomb()
pv1, pv2 = lat.prim_vecs
ysym = kwant.TranslationalSymmetry(pv2 - 2*pv1) # lattice symmetry in -y direction
xsym = kwant.TranslationalSymmetry(-pv2) # lattice symmetry in -x direction
# adatom lattice, for visualization only
adatom_lat = kwant.lattice.Monatomic(lat.prim_vecs, name='adatom')
def site_size(s):
return 0.25 if s.family in lat.sublattices else 0.35
def site_color(s):
return '#000000' if s.family in lat.sublattices else '#de962b'
def create_lead_h(W, symmetry, axis=(0, 0)):
lead = kwant.Builder(symmetry)
lead[lat.wire(axis, W)] = 0. * sigma_0
lead[lat.neighbors(1)] = kinetic
return lead
def create_lead_v(L,W, symmetry):
axis=(0, 0)
lead = kwant.Builder(symmetry)
lead[lat.wire(axis, W)] = 0. * sigma_0
lead[lat.neighbors(1)] = kinetic
return lead
def create_hall_cross(W, L, adatom_density=0.2, show_adatoms=False):
## scattering region ##
sys = kwant.Builder()
sys[lat.shape(cross(W, L), (0, 0))] = onsite_potential
sys[lat.neighbors(1)] = kinetic
adatoms = random_placement(sys, lat, adatom_density)
sys[(lat.hexagon(a) for a in adatoms)] = potential_shift
sys[(lat.hexagon_neighbors(a, 1) for a in adatoms)] = rashba
sys[(lat.hexagon_neighbors(a, 2) for a in adatoms)] = spin_orbit
if show_adatoms:
# no hoppings are added so these won't affect the dynamics; purely cosmetic
sys[(adatom_lat(*a) for a in adatoms)] = 0
## leads ##
leads = [create_lead_h(W, xsym)]
leads += [create_lead_v(L,W, ysym)]
leads += [lead.reversed() for lead in leads] # right and bottom leads
for lead in leads:
sys.attach_lead(lead)
return sys
def plot_ldos(sys, Es, params):
fsys = sys.finalized()
ldos = kwant.ldos(fsys, energy=Es, args=(params,))
ldos = ldos[::2] + ldos[1::2] # sum spins
kwant.plotter.map(fsys, ldos, vmax=0.1)
def plot_conductance(sys, energies,lead_i,lead_j, params):
fsys = sys.finalized()
data = []
for energy in energies:
smatrix = kwant.greens_function(fsys, energy, args=(params,))
data.append(smatrix.transmission(lead_i,lead_j))
plt.figure()
plt.plot(energies, data)
plt.xlabel("energy (t)")
plt.ylabel("conductance (e^2/h)")
plt.show()
def plot_spin_conductance(sys, energies,lead_i,lead_j, params):
fsys = sys.finalized()
data = []
for energy in energies:
smatrix = kwant.greens_function(fsys, energy, args=(params,))
ttdagger = smatrix._a_ttdagger_a_inv(lead_i, lead_j)
sigma_z_matrix = np.kron(np.eye(len(ttdagger)/2),sigma_z)
data.append(np.trace(sigma_z_matrix.dot(ttdagger)).real)
plt.figure()
plt.plot(energies, data)
plt.xlabel("energy (t)")
plt.ylabel("spin conductance (e^2/h)")
plt.show()
def solvematrix4(G,I0,I2,V1,V3):
# I=GxV
I1,I3,V0,V2=S('I1 I3 V0 V2'.split())
equations = [Eq(I0, G[0][0]*V0+G[0][1]*V1+G[0][2]*V2+G[0][3]*V3), \
Eq(I1, G[1][0]*V0+G[1][1]*V1+G[1][2]*V2+G[1][3]*V3), \
Eq(I2, G[2][0]*V0+G[2][1]*V1+G[2][2]*V2+G[2][3]*V3), \
Eq(I3, G[3][0]*V0+G[3][1]*V1+G[3][2]*V2+G[3][3]*V3)]
return solve(equations, [I1,I3,V0,V2], set=True, dict=False, rational=None, manual=True, minimal=True, particular=True, implicit=False, quick=True)
def plot_SH_angle(sys, energies, params):
fsys = sys.finalized()
data = []
I0=0
I2=0
V1=1.0
V3=0
for energy in energies:
smatrix = kwant.greens_function(fsys, energy, args=(params,))
G=np.zeros((4,4))
for i in range(4):
a=0
for j in range(4):
G[i,j]=smatrix.transmission(i, j)
if i != j: a+=G[i,j]
G[i,i]=-a
a=solvematrix4(G,I0,I2,V1,V3)[1]
result=list(a)
I1=result[0][0]
I3=result[0][1]
V0=result[0][2]
V2=result[0][3]
print(I1,I3,V0,V2)
Gspin=np.zeros((4,4))
for i in range(4):
a=0
for j in range(4):
ttdagger = smatrix._a_ttdagger_a_inv(i, j)
sigma_z_matrix = np.kron(np.eye(len(ttdagger)/2),sigma_z)
Gspin[i,j]=np.trace(sigma_z_matrix.dot(ttdagger)).real
if i != j: a+=Gspin[i,j]
Gspin[i,i]=-a
IS0=Gspin[1,0]*(V1-V0)+Gspin[2,0]*(V2-V0)+Gspin[3,0]*(V3-V0)
data.append(IS0/I1)
plt.figure()
plt.plot(energies, data)
plt.xlabel("energy (t)")
plt.ylabel("spin hall angle")
plt.show()
if __name__ == '__main__':
params = dict(gamma=1., ep=0, mu=0.1, V_I=0.007, V_R=0.0165) # Au adatoms
W=16
L=32
adatom_density=0.15
sys = create_hall_cross(W, L, adatom_density, show_adatoms=True)
kwant.plot(sys, site_color=site_color, site_size=site_size)
sys = create_hall_cross(W, L, adatom_density, show_adatoms=False)
plot_ldos(sys, 1e-5, params)
Es = np.linspace(-0.1, 0.1, 100)
plot_spin_conductance(sys, Es, 1, 0, params) # T01
plot_conductance(sys, Es, 1, 0, params) # TS01
plot_SH_angle(sys, Es, params)
