Exact Diagonalization of Fermi-Hubbard Models

exact-diagonalization-fermions

I’ll show you how to numerically diagonalize a fermionic lattice many-body quantum system. The technique has also been implemented in a python module called fermihubbard.py.

The Hamiltonian

We consider electrons that can move on a lattice. The electrons come in two flavors: spin up and spin down. Electrons are fermions and so they obey the Pauli exclusion principle, which means that no two fermions can reside in the same single lattice site.

Our system is described by the Fermi-Hubbard Hamiltonian,

\begin{align} \hat{H} &= \hat{H}_T + H_{\varepsilon} + \hat{H}_U \\ &= \sum_{\langle i, j \rangle} \sum_{\sigma \in \{\uparrow, \downarrow\}} \left( t_{ij} \, \hat{c}^\dagger_{\sigma i} \hat{c}^{\phantom{\dagger}}_{\sigma j} + \text{h.c.}\right) + \sum_{i} \varepsilon_i (\hat{n}_{\uparrow i} + \hat{n}_{\downarrow i}) + U \sum_{i} \hat{n}_{\uparrow i} \hat{n}_{\downarrow i}. \end{align}

It consists of three parts:

The hopping term, $\hat{H}_t$ describes the electrons moving between neighboring lattice sites $\langle i, j \rangle$ . This is described by the second quantization operators for creation, $\hat{c}^{\dagger}_{\sigma i}$ , and annihilation, $\hat{c}^{\phantom{\dagger}}_{\sigma j}$ , of an electron on site $i$ with spin $\sigma$ . Note that the hopping amplitudes must obey $t_{ij} = t^*_{ji}$ for hermiticity.
The electro-static term, $\hat{H}_\varepsilon$ , captures the effects of a varying electrostatic environment across different sites by adding a site-dependent energi, $\varepsilon_i$ . The operators are expressed through counting operators $\hat{n}_{\sigma i} = \hat{c}^{\dagger}_{\sigma i} \hat{c}^{\phantom{\dagger}}_{\sigma i}$ on the different sites, $i$ and spins ${\sigma = \{\uparrow,\,\downarrow\}}$ .
The Coulomb repulsion, $\hat{H}_U$ , describes the electrical repulsion (of magnitude $U$ ) between two electrons that reside on the same lattice site. Because the Pauli exclusion principle already prohibits two spin-up (or two spin-down) electrons to inhabit the same lattice site, this interaction is only present between spin-up and spin-down electrons.

It is the Coulomb repulsive term that complicates matters. If the interaction was absent, the electrons would not be able to “see” each other, and we could instead get away with just solving for the behavior of each electron separately. However, the interaction between the spin-up and spin-down electrons means we instead have to solve the full many-body quantum state.

The system

Let us label the lattices sites from $i = 1\ldots N$ . The electronic state can then be built from basis states representing a certain occupation of up-electrons and down-electrons. Consider as an example this basis state:

The fermionic many-body state: [u0, 00, ud, u0, 00, 0d, ..., ud, 0d] with "u" here references spin up and "d" spin down.

Formally it is equivalent to this second-quantization representation (notice the spin-up before spin-down convention to consistently account for the exchange statistics of fermions):

| {\uparrow} 0 - 00 - {\uparrow} {\downarrow} - 00 - 0{\downarrow} - \ldots - {\uparrow} {\downarrow} - 0 {\downarrow} \rangle = \hat{c}^{\dagger}_{\uparrow 0} \hat{c}^{\dagger}_{\uparrow 2} \hat{c}^{\dagger}_{\uparrow 3} \cdots \hat{c}^{\dagger}_{\uparrow (N-1)} \hat{c}^{\dagger}_{\downarrow 2} \hat{c}^{\dagger}_{\downarrow 5} \cdots \hat{c}^{\dagger}_{\downarrow (N-1)} \hat{c}^{\dagger}_{\downarrow N} | 0 \rangle.

The many-body states (like the ground state) that we are interested in, can then be represented as superpositions of many such basis states.

Simplification

The Hamiltonian does not contain terms that flip the spin of an electron (formally speaking it commutes with the counting operators $\hat{n}_{\uparrow}$ and $\hat{n}_{\downarrow}$ ). This means that we can separate the spin-up and the spin-down components of our state:

| {\uparrow} 0 - 00 - {\uparrow} {\downarrow} - 00 - 0{\downarrow} - \cdots - {\uparrow} {\downarrow} - 0 {\downarrow} \rangle = |{\uparrow} 0{\uparrow} 00 \cdots {\uparrow} 0 \rangle \otimes |00{\downarrow} 0{\downarrow} \cdots {\downarrow} {\downarrow} \rangle

Each of the $N_{\uparrow}$ spin-up electrons has $N$ sites to occupy and vice versa for the $N_{\downarrow}$ spin-down electrons. The total number of basis states is then given by a product of binomial coefficients:

n(N, N_{\uparrow}, N_{\downarrow})= \frac{N!}{(N-N_{\uparrow})!N_{\uparrow}!} \cdot \frac{N!}{(N-N_{\downarrow})!N_{\downarrow}!}.

This number grows quickly with system size, and means that an exact solution of the full many-body state is only feasible for small systems. Just to give you a sense of the scale here: For 20 sites with 20 electrons (10 up and 10 down) this amounts to the impressive

n(20, 10, 10) = (20!)^2 / (10!)^4 = 34.134.779.536.

Simply writing down a many-body state in this basis would amount to more than 34 billion complex numbers taking up more than 500 GB of data storage.

The Matrix

Even when limited to smaller systems, we must construct the Hamiltonian in matrix form. The standard technique is to insert two complete sets of basis states $|i\rangle$ and $|j\rangle$ , like so

\hat{H} = \sum_{i} |i \rangle \langle i | \: \hat{H} \: \sum_{j} | j \rangle \langle j | = \sum_{ij} \langle i | \hat{H} | j \rangle |i \rangle \langle j |.

To construct the matrix representation we need the amplitudes $\langle i | \hat{H} | j \rangle$ . The main challenge is the hopping, $\hat{H}_t$ , part of the Hamiltonian, which depends on the geometry of our lattice.

Consider e.g. the state $|{\uparrow}{\uparrow} 00 \rangle$ with two spin-up electrons on the first 2 of 4 total sites. The part of the hopping operator that moves an electron from site 0 to site 3 takes the form: $t_{03} \hat{c}^{\dagger}_{\uparrow 3} \hat{c}^{\phantom{\dagger}}_{\uparrow 0}$ . Applied to our example state, it would only give a non-zero amplitude only for the state where the electrons have actually moved, i.e. $\langle 0 {\uparrow} 0 {\uparrow} |$ . The actual amplitude is:

\begin{align} \langle 0 {\uparrow} 0 {\uparrow} | t_{03} \hat{c}^\dagger_{\uparrow 3} \hat{c}^{\phantom{\dagger}}_{\uparrow 0} |{\uparrow}{\uparrow} 00 \rangle &= t_{03} \langle 0 | \hat{c}^{\phantom{\dagger}}_{\uparrow 3} \hat{c}^{\phantom{\dagger}}_{\uparrow 1} \hat{c}^\dagger_{\uparrow 3} \hat{c}^{\phantom{\dagger}}_{\uparrow 0} \hat{c}^{\dagger}_{\uparrow 0} \hat{c}^{\dagger}_{\uparrow 1} | 0 \rangle \\ &= - t_{03} \end{align}

where the minus sign arises due to the fermionic exchange statistics.

Lin tables

In order to construct the matrix we need to number the basis states. Usually this is handled by Lin tables (a hash table). The idea is to translate each spin-component using bit-patterns, so they can be represented with binary numbers, so,

|{\uparrow}{\uparrow} 00 \rangle \sim 1100 \ \rightarrow 12.

With 4 sites and 2 $\uparrow$ electrons, there are a total of 6 $\uparrow$ corresponding with this Lin table.

State number	State	Hash
0	0011	3
1	0101	5
2	0110	6
3	1001	9
4	1010	10
5	1100	12

This table allows us to translate freely between the state (i.e. 0110) and the state number (i.e. 2).

Square Example

Consider, as an example, four sites connected in a ring-geometry with equal hoppings.

The 4-site ring example

We are often interested in the half-filled zero-spin case, which in this example correspond to a filling with $2 \uparrow$ and $2 \downarrow$ electrons. For the $2\uparrow$ (or $2 \downarrow$ ) sub-space the resulting hopping Hamiltonian becomes:

\hat{H}_{t, 2\uparrow} = t \left( \begin{array}{cccccc} 0 & -1 & 0 & 0 & 1 & 0 \\ -1 & 0 & -1 & -1 & 0 & 1 \\ 0 & -1 & 0 & 0 & -1 & 0 \\ 0 & -1 & 0 & 0 & -1 & 0 \\ 1 & 0 & -1 & -1 & 0 & -1 \\ 0 & 1 & 0 & 0 & -1 & 0 \\ \end{array} \right)

The full hopping Hamiltonian, can then be constructed as the Kronecker product:

\hat{H}_{t} = \hat{H}_{t, 2\uparrow} \otimes \hat{H}_{t, 2 \downarrow}

The interaction Hamiltonian, $\hat{H}_U$ , is diagonal and can be computed directly by applying logical AND between the spin-up and spin-down bit patterns. Similarly the electro-static term, $\hat{H}_{\varepsilon}$ , is computed as the sum of the spin-up and spin-down bit patterns.

The fermihubbard module

I have implemented an algorithm for constructing the Hamiltonian in a python module called fermihubbard.py. It requires the single particle representation of the hopping and electrostatic terms, and a similar matrix describing the Coulomb term.

For our 4-ring example:

1
import numpy
2
import fermihubbard
3

4
# The onsite energies (all zero)
5
H1E = np.zeros((4, 4))
6

7
# Hopping on a ring (amplitude t = 1)
8
H1T = np.zeros((4, 4))
9
for i in xrange(4):
10
    H1T[i, (i+1) % 4] = 1.
11
    H1T[(i+1) % 4, i] = 1.
12

13
# Interaction (here onsite with U = 2)
14
H1U = 2 * numpy.eye(4)
15

16
# Construct the model
17
model = fermihubbard.Model(H1E, H1T, H1U)

We can restrict ourselves to the relevant subspace:

1
# Construct a specific chargestate with eight electrons:
2
number_sector_4 = m.numbersector(4)
3

4
# Consider the sector with net zero spin in the z-direction, meaning that
5
# 4 electrons = 2 spin-up electrons + 2 spin-down electrons
6
spin_component_0 = number_sector_4.szsector(0)

The Lin table and Hamiltonian can be generated for the relevant subspace:

1
# Let us then take a look at the basis. Generate the Lin-table:
2
lin_table = lin.Table(4, 2, 2)
3

4
# Compute the Hamiltonian terms
5
(HTu, HTd, HUd) = spin_component_0.hamiltonian

From the three parts, the complete Hamiltonian can be constucted in e.g. sparse matrix format that is also suited for larger systems.

1
import scipy.sparse as sparse
2

3
# The total Hamiltonian can be generated from these parts.
4
Nu, Nd = lin_table.Ns
5

6
# The interaction part in sparse format
7
HU = sparse.coo_matrix(
8
    (HUd, (np.arange(Nu * Nd), np.arange(Nu * Nd))),
9
    shape=(Nu * Nd, Nu * Nd)
10
).tocsr()
11

12
# The kronecker product of the two hopping sectors.
13
HT = (
14
    sparse.kron(np.eye(Nd), HTu, 'dok') +
15
    sparse.kron(HTd, np.eye(Nu), 'dok')
16
)
17

18
# The total Hamiltonian
19
H = HU + HT

Diagonalization

For small systems the Hamiltonian can be diagonalized fully using standard algorithms. As the system grows the iterative Lanczos method can bed used to only calculate the low energy states.

1
# The 10 lowest energy states using an iterative eigensolver.
2
from scipy.linalg import eigh
3
energies, states = eigh(H, subset_by_index=(0, 10))
4

5
print('Lowest eigenvalue: ', energies[0])

Conclusion

To sum up: We have “solved” fermi-Hubbard models using exact diagonalization. The resulting energies and eigenstates can then be used to compute various observables.

An advanced example is e.g. the cotunneling current in molecular junstions as done in this paper from my PhD.