Topological Hall effect in four terminal devices

KWANT script below describes a skyrmion and computes the Hall Resistance and Spin Hall angle of the four-terminal device.

from types import SimpleNamespace
from matplotlib import pyplot
from math import cos, sin, pi
import numpy as np
import scipy.stats as reg
import kwant
from sympy import S, Eq, solve, Symbol

lat = kwant.lattice.square()

s_0 = np.identity(2)
s_z = np.array([[1, 0], [0, -1]])
s_x = np.array([[0, 1], [1, 0]])
s_y = np.array([[0, -1j], [1j, 0]])

def HedgeHog(site,ps):
        x,y = site.pos
        r = ( x**2 + y**2 )**0.5
        theta = (np.pi/2)*(np.tanh((ps.r0 - r)/ps.delta) + 1)
        if r != 0:
            Ex = (x/r)*np.sin(theta)*s_x + (y/r)*np.sin(theta)*s_y + np.cos(theta)*s_z
        else:
            Ex = s_z
        return 4*s_0 + ps.Ex * Ex
       

def Lead_Pot(site,ps):
    return  4*s_0 +  ps.Ex * s_z

def MakeSystem(ps, show = False):
    H  = kwant.Builder()

    def shape_2DEG(pos):
        x,y = pos
        return  ( (abs(x) < ps.L) and (abs(y) < ps.W) ) or ( 
            (abs(x) < ps.W) and (abs(y) < ps.L))
    
    H[lat.shape(shape_2DEG,(0,0))] = HedgeHog
    H[lat.neighbors()] = -s_0
    
    # Leads
    sym_x = kwant.TranslationalSymmetry((-1,0))
    H_lead_x = kwant.Builder(sym_x)
    shape_x = lambda pos: abs(pos[1])<ps.W and pos[0]==0 
    H_lead_x[lat.shape(shape_x,(0,0))] = Lead_Pot
    H_lead_x[lat.neighbors()] = -s_0
   
    sym_y = kwant.TranslationalSymmetry((0,-1))
    H_lead_y = kwant.Builder(sym_y)
    shape_y = lambda pos: abs(pos[0])<ps.W and pos[1]==0 
    H_lead_y[lat.shape(shape_y,(0,0))] = Lead_Pot 
    H_lead_y[lat.neighbors()] = -s_0
    
    H.attach_lead(H_lead_x)
    H.attach_lead(H_lead_y)
    H.attach_lead(H_lead_y.reversed())
    H.attach_lead(H_lead_x.reversed())
    
    if show:
        kwant.plot(H)

    return H
 

def Transport(Hf,EE,ps):
    smatrix = kwant.smatrix(Hf, energy=EE, args=[ps])
    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 
        
    V = np.linalg.solve(G[:3,:3], [1.,0,0])
    Hall = V[2] - V[1]
    
    return G, Hall
	
ps = SimpleNamespace(L=45, W=40, delta=10, r0=20, Ex=1.)

H = MakeSystem(ps, show=True)
Hf = H.finalized()

def Vz(site):
    Hd = HedgeHog(site,ps)
    return (Hd[0,0] - Hd[1,1]).real 

def Vy(site):
    Hd = HedgeHog(site, ps)
    return Hd[0,1].imag 

kwant.plotter.map(H, Vz);
kwant.plotter.map(H, Vy);

# Hall Resistance
ps = SimpleNamespace(L=20, W=15, delta=3, r0=6, Ex=1.)

H = MakeSystem(ps, show=False)
Es = np.linspace(0.1,5.,50)
Hf = H.finalized()
dataG , dataHall = [],[]

for EE in Es: 
    ps.delta = EE
    energy = 2.
    G,Hall = Transport(Hf, energy, ps)
    dataHall.append(Hall)

pyplot.plot(Es, dataHall, 'o-', label="Skyrmion")
pyplot.xlabel('Domain width $\delta$')
pyplot.ylabel('Hall Resistance')
pyplot.title('Topologial Hall  Resistance')
pyplot.legend();
pyplot.show()

# Spin Hall angle
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(fsys, energies, ps):
    data = []

    I0=0
    I2=0
    V1=1.0
    V3=0

    for energy in energies:
        smatrix = kwant.greens_function(fsys, energy, args=[ps])
        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),s_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)

    pyplot.figure()
    pyplot.plot(energies, data)
    pyplot.xlabel("energy (t)")
    pyplot.ylabel("spin hall angle")
    pyplot.show()
plot_SH_angle(Hf, Es, ps)