Simulation of the 1D Tight-Binding Model, full tutorial https://mareknarozniak.com/2020/05/07/tight-binding/

run.py

import matplotlib.pyplot as plt
import numpy as np
from qutip import basis


def _ket(n, N):
    return basis(N, n)


def _bra(n, N):
    return basis(N, n).dag()


# spacing between atoms
a = 2

# number of atoms
N = 23

# total length of the 1D chain
L = N*a

# onsite energy
eps = 0.25

# hopping matrix element
t = 0.5

# construct the Hamiltonian
H  = sum([eps*_ket(n, L)*_bra(n, L) for n in range(0, L)])
H -= sum([t*_ket(n, L)*_bra(n + 1, L) for n in range(0, L - 1)])
H -= sum([t*_ket(n, L)*_bra(n - 1, L) for n in range(1, L)])

if N <= 8:
    print('Hamiltonian')
    print(H)

evals, ekets = H.eigenstates()

if N <= 8:
    print()
    print('Eigenvalues')
    print(evals)

# satisfy periodic boundary conditions we need E(-k) = E(k)
numerical = np.concatenate((np.flip(evals, 0), evals), axis=0)

k = np.linspace(-np.pi/a, np.pi/a, 2*L)
exact = eps - 2.*t*np.cos(k*a)

xticks = np.linspace(-np.pi/a, np.pi/a, 9)
xlabels = ['' for k in xticks]
xlabels[0] = '$-\\frac{\pi}{a}$'
xlabels[-1] = '$\\frac{\pi}{a}$'

fig, axs = plt.subplots()
axs.set_xlim(-np.pi/a, np.pi/a)
axs.set_title('Tight-Binding Model, $N='+str(N)+', a='+str(a)+'$')
axs.set_ylabel('$E$')
axs.set_xlabel('$ka$')
axs.axvline(x=0., color='k')
axs.plot(k, numerical, 'ro', label='Eigenvalues of H')
axs.plot(k, exact, label='$E(k) = \epsilon_0 - 2 t \cos(ka)$')
axs.set_yticks([eps-2.*t, eps-t, eps, eps+t, eps+2.*t])
axs.set_yticklabels(['$\epsilon_0-2t$', '$\epsilon_0-t$', '$\epsilon_0$', '$\epsilon_0+t$', '$\epsilon_0+2t$'])
axs.set_xticks(xticks)
axs.set_xticklabels(xlabels)
axs.legend()
axs.grid(True)

# fig.savefig('tight_binding_theory.png')
fig.savefig('tight_binding_diag.png')

tight_binding_diag.png

tight_binding_theory.png