This post is available as an IPython notebook.
Hamiltonians on one dimensional chains provide a good sanity check when trying out new numerical methods. The Jordan-Wigner transform is an easy way to obtain the eigenvalues of systems of fermions. For an annihiliation operator 
, the transform is defined as 
, where 
is the Pauli operator along the 
axis, and 
is the operator 
acting on site 
.
The series of Kronecker products would be an ideal application of a nest function, but Python does not have such a function built in. Instead, we define a recursive function to get the desired Kronecker product:
import numpy as np
def nested_kronecker_product(a):
if len(a) == 2:
return np.kron(a[0],a[1])
else:
return np.kron(a[0], nested_kronecker_product(a[1:]))
With this function, the transform is easy to define. We pad the transform with identity operators to match the lattice length:
def jordan_wigner_transform(j, lattice_length):
sigma = np.array([[0, 1], [0, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
I = np.eye(2)
operators = []
for k in range(j):
operators.append(sigma_z)
operators.append(sigma)
for k in range(lattice_length-j-1):
operators.append(I)
return -nested_kronecker_product(operators)
Now we can get our transformed fermionic operators:
lattice_length = 4
a = []
for i in range(lattice_length):
a.append(jordan_wigner_transform(i, lattice_length))
Next we define a Hamiltonian of interest, namely spinless fermions on an open chain:
where 
and 
indicate neighbors on the chain. The Python definition of this Hamiltonian is
def hamiltonian(gam, lam, a, lattice_length):
H = 0
for i in range(lattice_length - 1):
H += a[i].T.dot(a[i+1]) - a[i].dot(a[i+1].T)
H -= gam*(a[i].T.dot(a[i+1].T) - a[i].dot(a[i+1]))
for i in range(lattice_length):
H -= 2*lam*(a[i].dot(a[i].T))
return H
Here we observed the fermionic canonical commutation relations and performed substitutions. To get the eigenvalues for a particular parameter combination, we write
gam, lam =1, 1
H = hamiltonian(gam, lam, a, lattice_length)
eigenvalues = np.linalg.eig(H)[0]
sorted(eigenvalues)
