Dear Sir,

I am trying to plot the energy as a function of magnetic field for a 3D case, but I am getting tedious results. The system is infinite in two directions and has some width in the third direction. Please have a look at the code attached below. I tried a lot but failed. Is the code correct or I am wrong somewhere.

Thank you.

import kwant
import scipy.sparse.linalg as sla
import matplotlib.pyplot as plt
import tinyarray
import numpy as np
from numpy import cos, sin, pi
import cmath
from cmath import exp

sigma_0 = tinyarray.array([[1, 0], [0, 1]])
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])

def make_system(a=1, L=30, W=10, H=10, t=1.0, t_x=1.0, t_y=1.0, t_z=1.0, lamda=0.1, beta=1.05):
def onsite(site):
return (t_z * cos(beta) + 2 * t) * sigma_z

def hoppingx(site0, site1):
return (-0.5 * t * sigma_z - 0.5 * 1j * t_x * sigma_x)

def hoppingy(site0, site1):
return -0.5 * t * sigma_z - 0.5 * 1j * t_y * sigma_y

def hoppingz(site0, site1, B):
y = site1.pos[1]
return (-0.5 * t_z * sigma_z - 0.5 * 1j * lamda * sigma_0) * exp(2 * pi * 1j * B * a * (y-40))

syst = kwant.Builder()
lat = kwant.lattice.cubic(a)
syst[(lat(z, y, x) for z in range(H) for y in range(W) for x in range(L))] = onsite
syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx

lead1=kwant.Builder(kwant.TranslationalSymmetry((0,-a,0)))
lead1[(lat(z,y,x) for z in range(H)for y in range(W)for x in range(L))]=onsite
lead1[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
lead1[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
lead1[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx

syst.attach_lead(lead1)
syst.attach_lead(lead1.reversed())

lead2=kwant.Builder(kwant.TranslationalSymmetry((-a,0,0)))
lead2[(lat(z,y,x) for z in range(H)for y in range(W)for x in range(L))]=onsite
lead2[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
lead2[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
lead2[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx

syst.attach_lead(lead2)
syst.attach_lead(lead2.reversed())
syst = syst.finalized()
return syst

def analyze_system():
syst = make_system()
kwant.plot(syst)
Bfields = np.linspace(0, 0.002, 100)
energies = []
for B in Bfields:
ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
ev, evec = sla.eigsh(ham_mat.tocsc(), k=20, sigma=0)
energies.append(ev)
#print(energies)

plt.figure()
plt.plot(Bfields, energies)
plt.xlabel("magnetic field [${10^-3 h/e}$]")
plt.ylabel("energy [t]")

plt.ylim(0, 0.11)
plt.show
def main():
syst = make_system()
analyze_system()
main()

With Best Regards
NAVEEN YADAV
Ph.D Research Scholar
Deptt. Of Physics & Astrophysics
University Of Delhi.