# Quantum statistics and Hamiltonians

Let us first note that the Hamiltonian in the Schrödinger equation

$H \Psi = E \Psi$

is independent of the symmetry of the wavefunction. For a many-particle system, the distinguishable case makes no symmetry assumptions on Ψ, and it follows that the spectrum of the distinguishable-particle system contains the spectrum of the symmetric (bosonic) and the anti-symmetric (fermionic) system. It contains also states of mixed symmetry, neither bosonic nor fermionic. In an earlier lecture, we saw how the symmetry requirement was introduced into the density matrix via a sum over permutations. Here, we see how the symmetry requirement is realized in first and second quantization for lattice systems of distinguishable particles, of bosons, and of fermions.

NB: In this lecture, particles (even fermions) are spinless.

## An exact diagonalization algorithm for distinguishable quantum lattice particles

Here is an exact diagonalization algorithm for the one-dimensional Hubbard model of N spin-less distinguishable particles on L lattice sites with periodic boundary conditions, with a potential energy
$\sum_k U n_k(n_k-1)$
and a hopping term -t, if there is one particle hopping from one site to one of its neighbors.
import numpy
def general_count(i_min,i_max,N):
#
# enumerate states: each particle can sit on any of the N sites of the lattice
#
I_vec = [i_min] * N
states=[]
condition = False
while condition == False:
states.append(I_vec[:])
for k in range(N-1,-2,-1):
if k == -1: condition = True
if I_vec[k] == i_max:
I_vec[k] = i_min
else:
I_vec[k] += 1
break
return states

U = 4
t = 0.75
L = 4
N = 3
n_states = L**N
a_mat = numpy.zeros((n_states,n_states))
nbr = {}

nbr[0] = [L-1,1]
nbr[L-1] = [L-2,0]
for i in range(1, L - 1):
nbr[i] = (i - 1, i + 1)
print nbr
states = general_count(0, L - 1, N)
print states
for i in states:
ii = states.index(i)
a_mat[ii][ii] = U * sum(i.count(k) * (i.count(k) - 1) for k in range(L))
for k in range(len(i)):
for l in nbr[i[k]]:
j= i[:]
j[k] = l
print i, j
jj = states.index(j)
a_mat[ii][jj] = -t
print a_mat
w, v = numpy.linalg.eigh(a_mat)
print " Spectrum of H"
for x in w:
print x

Here is the (partial) list of the 64 states generated by the above program:
[0, 0, 0]
[0, 0, 1]
[0, 0, 2]
[0, 0, 3]
[0, 1, 0]
[0, 1, 1]
[0, 1, 2]
...

Here is the (partial) 64 × 64 Hamiltonian matrix:
[[ 24.    -0.75  0.   ..., 0. 0. 0. ]
[  -0.75   8.   -0.75 ..., 0. 0. 0. ]
[   0.    -0.75  8.   ..., 0. 0. 0. ]
...,
[ 0. 0. 0. ..., 8. -0.75 0. ]
[ 0. 0. 0. ..., -0.75 8. -0.75]
[ 0. 0. 0. ..., 0. -0.75 24. ]]

Here, the element H00 = 24, because three particles are on one site (0), and
$H_{01} = -3/4$
because particle no 3 hops from site 0 to site 1.

## An exact diagonalization algorithm for quantum lattice bosons

Here is an equivalent exact diagonalization algorithm for the one-dimensional Hubbard model of spin-less bosonic particles. Second quantized
notations are used, so that the hamiltonian is now

$H = - t \sum_{\langle i,j\rangle}c_i^+ c_j + \text{h.c.} + U \sum_i n_i (n_i-1)$
We note the combinatorial factors
$c_i^+ |n_i\rangle = \sqrt{n_i+1} |(n_i+1)\rangle$
and
$c_i |n_i\rangle = \sqrt{n_i} |(n_i-1)\rangle$
(see R. P. Feynman ''Statistical mechanics: a set of lectures'' (Benjamin Cummings, 1972) to understand the square roots). Besides the diagonal terms, the Hamiltonian connects only states with correspond to a single hop of a particle from site i to a neighboring site j.

import numpy
from math import sqrt
def general_count(i_min, i_max, L,N):
I_vec = [i_min] * L
states = []
condition = False
while condition == False:
if sum(I_vec) == N: states.append(I_vec[ : ])
for k in range(L-1, -2, -1):
if k == -1: condition = True
if I_vec[k] == i_max:
I_vec[k] = i_min
else:
I_vec[k] += 1
break
return states
U = 4.0
t = 0.75
L = 4
N = 3
nbr = {}
nbr[0] = [L - 1, 1]
nbr[L - 1] = [L - 2, 0]
for i in range(1, L - 1):
nbr[i] = (i-1, i + 1)
print nbr
states = general_count(0, N, L, N)
a_mat = numpy.zeros((len(states), len(states)))
for i in states:
ii = states.index(i)
a_mat[ii][ii] = U * sum(k * (k-1) for k in i)
for k in range(L):
for l in nbr[k]:
j = i[ : ]
if (j[l] > 0):
hopping = sqrt((i[k] + 1) * i[l])
j[k] = i[k] + 1
j[l] = i[l] - 1
print j
jj = states.index(j)
a_mat[ii][jj] = -hopping * t
print a_mat
w, v = numpy.linalg.eigh(a_mat)
for x in w:
print x

Here, the states have length 4, and there are 20 of them ( see W. Krauth ''Statistical mechanics: Algorithms and Computations'' (Oxford University Press, 2006) See also [http://www.smac.lps.ens.fr/ the book's website] eq. 4.12)
[0, 0, 0, 3]
[0, 0, 1, 2]
[0, 0, 2, 1]
[0, 0, 3, 0]
[0, 1, 0, 2]
[0, 1, 1, 1]
[0, 1, 2, 0]
...


and the matrix is
[ 24.        -1.29903811 0.0         0.0        0.0        0.0 ...
[ -1.29903811 8.0       -1.5         0.0       -0.75       0.0 ...
[  0.0       -1.5        8.0        -1.29903811 0.0       -1.06066017 ...
[  0.0        0.0       -1.29903811 24.0        0.0        0.0 ...
[  0.0       -0.75       0.0         0.0        8.0       -1.06066017 ...
[  0.0        0.0       -1.06066017  0.0       -1.06066017 0.0 ...
...

here, we note that
$1.299038.. = \frac{3}{4} \sqrt{3}$

## Exact diagonalization algorithm for (spinless) lattice fermions

The Hubbard hamiltonian for lattice fermions is usually written for two spin species, and
For fermions, the program is much simpler, because there are only 4 states (for 3 particles on 4 sites). Here is the program computing the Hamiltonian and the spectrum:

import numpy
import numpy.linalg
U = 4.0
t = 0.75
L = 4
n_states = 4
a_mat = [[0., -t, 0., -t], [-t, 0., -t, 0.], [0., -t, 0., -t], [-t, 0., -t, 0.]]
print a_mat
w, v = numpy.linalg.eigh(a_mat)
for dummy in w:
print dummy

Fermionic states are
(1, 1, 1, 0)
(1, 1, 0, 1)
(1, 0, 1, 1)
(0, 1, 1, 1)

The Hamiltonian matrix is
[[ 0.0, -0.75, 0.0, -0.75],
[ -0.75, 0.0, -1.0,  0.0],
[  0.0, -0.75, 0.0, -0.75],
[ -0.75, 0.0, -0.75, 0.0]]

The eigenvalues of this matrix are -1.5, 0, 0, 1.5

## Comparing spectra

600px|left|visualizing two spectra with vimdiff Here we visualize the spectra of hubbard_bose.py (to the left) and of hubbard_dist.py (to the right), using an editor (vimdiff) capable of comparing two files. We see that the groundstate is of bosonic symmetry, but that not all of the states of the distinguishable hamiltonian appear either for the bosons or for the fermions.

## Checking states

We discussed before that the groundstate of the distinguishable-particle problem had bosonic statistics, so let us check that this is indeed true:
-2.28331340813 Energy... (this is the distinguishable-particle
and the bosonic groundstate)
(0.0067261131197429158, [1, 1, 1])
(0.0067261131197430841, [3, 3, 3])
(0.0067261131197431639, [0, 0, 0])
(0.006726113119743274, [2, 2, 2])
cc
(0.039285453143280247, [0, 1, 0])
(0.039285453143280691, [0, 0, 3])
(0.03928545314328085, [3, 0, 0])
(0.039285453143280878, [1, 0, 0])
(0.039285453143280892, [0, 3, 0])
(0.03928545314328092, [1, 1, 0])
(0.03928545314328092, [3, 0, 3])
(0.039285453143280927, [0, 3, 3])
(0.039285453143280927, [3, 2, 3])
(0.039285453143280968, [1, 1, 2])
(0.039285453143280975, [2, 1, 1])
(0.039285453143280982, [2, 3, 3])
(0.039285453143280989, [1, 0, 1])
(0.039285453143280989, [3, 3, 0])
(0.039285453143280989, [3, 3, 2])
(0.039285453143280996, [1, 2, 1])
(0.039285453143281024, [0, 1, 1])
(0.039285453143281059, [2, 2, 3])
(0.039285453143281086, [2, 1, 2])
(0.039285453143281086, [2, 2, 1])
(0.0392854531432811, [1, 2, 2])
(0.039285453143281107, [3, 2, 2])
(0.039285453143281128, [2, 3, 2])
(0.039285453143281177, [0, 0, 1])
...etc

and let's check, for fun, the fermionic statistics of the state with energy -1.5:
==
-1.5 (this is the fermionic groundstate, and an excited state of the distinguishable-particle Hamiltonian)
==
(-0.20412414523193678, [3, 0, 2])
(-0.20412414523193539, [1, 2, 3])
(-0.20412414523193423, [3, 1, 2])
(-0.20412414523193395, [1, 2, 0])
(-0.20412414523193284, [1, 3, 0])
(-0.20412414523193179, [3, 0, 1])
(-0.20412414523193165, [0, 1, 2])
(-0.20412414523192998, [0, 2, 3])
(-0.20412414523192959, [0, 1, 3])
(-0.20412414523192821, [2, 0, 1])
(-0.20412414523192746, [2, 3, 0])
(-0.20412414523192479, [2, 3, 1])
(0.20412414523191447, [1, 0, 2])
(0.20412414523191599, [1, 0, 3])
(0.20412414523192099, [1, 3, 2])
(0.20412414523192221, [2, 0, 3])
(0.20412414523193054, [2, 1, 3])
(0.2041241452319319, [0, 3, 2])
(0.20412414523193434, [2, 1, 0])
(0.20412414523193689, [3, 1, 0])
(0.20412414523193881, [3, 2, 0])
(0.20412414523194106, [0, 3, 1])
(0.2041241452319435, [3, 2, 1])
(0.20412414523194503, [0, 2, 1])

## Conclusion

Specifically, we have studied the Hubbard model (for 3 particles on 4 sites with periodic boundary conditions)
1. for ''distinguishable particles'' (64 × 64 Hamiltonian matrix, first program)
2. for ''bosons'' (20 × 20 Hamiltonian using creation and annihilation operators, second program)
3. for ''fermions'' (4 × 4 Hamiltonian matrix, final program).

The most important observation was that the spectrum of the distinguishable-particle Hamiltonian contained all possible states, the symmetric ones, the antisymmetric ones, but also the states with mixed symmetry. We also notice that the ground state of the Fermionic problem, in fact, is an excited state of the distinguishable-particle Hamiltonian. This, in a sense, is the origin of the '''fermion sign problem'''. It is interesting to see how three different formulations yield different aspects of the same system: the quantum many-body problem. We restricted ourselves to a very small number of sites and of particles, because of the exponential increase in the size of the Hamiltonian matrix.

More specifically, we have performed a first socalled '''exact diagonalization''' of a quantum hamiltonian. Exact diagonalization is a mature field of research which allows to compute the spectrum, or parts of the spectrum (for example the ground state only) of lattice hamiltonians. For thermodynamics, it is mostly used for fermions. Bosons can be handled very successfully using Monte Carlo algorithms, so that exact diagonalization methods are less common. The main problem in exact diagonalization is the '''exponential growth''' of the Hilbert space, which seriously limits the system sizes which can be taken into account. The main advantage is, as the name indicates, that results are numerically exact, without error bars, so that '''extrapolation''' can be very powerful.

# Lattice Monte Carlo algorithms

In earlier lectures, we treated the case of Monte Carlo algorithms in the continuum. Key elements were the convolution property, the explicit solution for the free density matrix, and the Trotter formula. These elements remain valid for the Monte Carlo

## Lattice bosons

For lattice bosons, the procedure is the same as for the continuum system,

## Lattice Fermions

For lattice fermions, we will treat a successful algorithm next time...

## Quantum spins

... unfortunately not treated in this course...Quantum Monte Carlo, similar to what we learned earlier, is possible...

### Computation of Fermion determinant

The formula of this section can be derived using Grassmann algebras. We prefer a pedestrian derivation, provided by J. E. Hirsch in the appendix of his famous paper J. E. Hirsch ''Two-dimensional Hubbard model: Numerical simulation study'' Physical Review B 31, 4403

$Z = \sum_{\text{states} \alpha } \langle \alpha |\exp(- \beta H) | \alpha \rangle$
We first show how to treat bilinear terms, such as

$\text{Tr} \exp( - c_i^+ B_{ij} c_j) = \text{Tr} \prod_\mu \exp( - c_\mu^+ b_\mu c_\mu)$
For this simple case, we can expand the exponential into

$\exp( - c_\mu^+ b_\mu c_\mu) = 1 - c_\mu^+ b_\mu c_\mu + \frac{1}{2}c_\mu^+ b_\mu c_\mu c_\mu^+ b_\mu c_\mu - \frac{1}{3!}.....$
Using c_\mu^+ c_\mu + c_\mu c_\mu^+ = 1, we reach

$\exp( - c_\mu^+ b_\mu c_\mu) = 1 [[[+ c_\mu^+ b_\mu c_\mu ]]] + \frac{1}{2}c_\mu^+ b_\mu^2 c_\mu - \frac{1}{3!}c_\mu^+ b_\mu^3 c_\mu .... [[[-c_\mu^+ b_\mu c_\mu ]]] = \prod_\mu [ 1 + ( \exp (- b_\mu) - 1) c_\mu^+ c_\mu$
This expression is easily evaluated. It yields the expression
$\det[1 + \exp(-B)]$
More generally, we find the expression:

$\text{Tr} \exp( - c_i^+ A_{ij} c_j)\exp( - c_i^+ B_{ij} c_j) = \det[1 + \exp(-A) \exp(-B) ]$

which can be derived from the action on single-particle states, or via Grassmann algebras. This non-trivial formula is the key to determinantal fermion methods.