I'm confused by some problem when doing time-dependent spin related transport calculations by tkwant. Let me briefly introduce what the main problem is. I'm now considering a 1D chain model with some magnetic spin moment on some of the sites. By applying a
voltage to one lead, I want to know how the spin accumulate near the other lead. My problem is that in my code the calculated spin density and charge density near the second lead are always same which I believe should not be true because this means that only
one spin channel is always occupied and the other is always empty! My code is as follows:
import tkwant
import kwant
from matplotlib import pyplot as plt
from math import pi, sqrt, cos, sin
from cmath import atan
import numpy as np
############ spin matrix #############
s0=np.array([[1,0],[0,1]])
sx=np.array([[0,1],[1,0]])
sy=np.array([[0,-1j],[1j,0]])
sz=np.array([[1,0],[0,-1]])
length=10
l_M=1
M=0.2
t=1
times=range(200)
lat = kwant.lattice.general([(5,0,0),(0,5,0),(0,0,1)],
[(2.5,2.5,0)],norbs=2)
def scattering_region(pos):
x,y,z=pos
return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z>=0 and z<=length-1
def mag_region0(pos):
x,y,z=pos
return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z<=-1 and z>=-l_M
def mag_region1(pos):
x,y,z=pos
return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z>=length and z<=length+l_M-1
syst=kwant.Builder()
syst[lat.shape(scattering_region, (2.5, 2.5, 0))]=0*s0
syst[lat.shape(mag_region0, (2.5, 2.5, -1))]=M*sz
syst[lat.shape(mag_region1, (2.5, 2.5, length))]=M*sz
for i in range (length+2*l_M-1):
syst[lat.sublattices[0](0,0,i-l_M),lat.sublattices[0](0,0,i-l_M+1)]=t*s0
def lead_region(pos):
x,y,z=pos
return abs(x-2.5)<=1 and abs(y-2.5)<=1
lead=kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0,0,-1))),conservation_law=-sz)
lead[lat.shape(lead_region, (2.5,2.5,0))]=0*s0
lead[lat.neighbors()]=t*s0
def phi(time):
Vb, tau= 0.6, 20
if time > tau:
return Vb*time - Vb*tau/2
return Vb*time*time/2/tau
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
kwant.plot(syst,site_size=0.2,site_lw=0.01,hop_lw=0.05)
tkwant.leads.add_voltage(syst,0,phi)
syst=syst.finalized()
def current_cut(site_to, site_from):
return site_from.pos[2]<=length-1 and site_to.pos[2]>length-1
def density_cut(site):
return site.pos[2]>=length and site.pos[2]<=length+l_M-1
current_operator=kwant.operator.Current(syst,where=current_cut)
charge_density_operator=kwant.operator.Density(syst,s0,where=density_cut)
spin_density_operator=kwant.operator.Density(syst,sz,where=density_cut)
occupy=tkwant.manybody.lead_occupation(chemical_potential=-0.5)
print('make up of state ...')
state=tkwant.manybody.State(syst,tmax=max(times),occupations=occupy,refine=False,combine=True)
Currents=[]
Charge=[]
Spin=[]
for time in times:
print('%.4f percent ...' % (100*times.index(time)/len(times)))
state.evolve(time)
current=state.evaluate(current_operator)
charge=state.evaluate(charge_density_operator)
spin=state.evaluate(spin_density_operator)
Currents.append(current)
Charge.append(charge)
Spin.append(spin)
np.savetxt('Times_linear_chain.txt',times)
np.savetxt('Current_linear_chain.txt',Currents)
np.savetxt('Charge_density_linear_chain.txt',Charge)
np.savetxt('Spin_density_linear_chain.txt',Spin)
plt.figure()
plt.plot(times,Currents)
plt.xlabel('time $t$')
plt.ylabel('Current $I$')
plt.show()
plt.figure()
plt.subplot(211)
plt.plot(times,Charge)
plt.xlabel('time $t$')
plt.ylabel('Charge density')
plt.subplot(212)
plt.plot(times,Spin)
plt.xlabel('time $t$')
plt.ylabel('Spin density')
plt.show()