# Is there a simple way to numerically solve interacting fermion systems like for spins?

I know that if one has a spin chain of N spins, with the interaction described by some Hamiltonian, one can find the eigenstates and energies by using Kronecker products to write the Hamiltonian as a $$(2s+1)^N \times (2s+1)^N$$ matrix and then digonalizing that matrix.

For example, the 1D quantum Ising model:

$$H = - \sum_i(J\sigma_i^x\sigma_{i+1}^x+g\sigma_i^z)$$

Written out explicitly, where each term is a $$2^N \times 2^N$$ matrix:

$$H = -J(\sigma^x \otimes \sigma^x \otimes I \otimes\cdots \otimes I+ I \otimes \sigma^x \otimes \sigma^x \otimes I \otimes\cdots \otimes I+ \cdots) - g(\sigma^z \otimes I \otimes \cdots \otimes I + \cdots)$$

Of course as the system size increases this matrix quickly becomes impossible to diagonalize, but at least the idea is clear.

Does there exist an analogous simple procedure for exactly solving a fermion system given a second quantized Hamiltonian?

I've seem examples of how to solve e.g. the Hubbard model numerically, but it seems like one needs a very ad hoc sort of procedure.

(Addendum: ok, I see how one can solve a fermion system pretty much the same way by having particle/no particle $$\leftrightarrow$$ spin up/spin down. But this only works if particle number is not conserved? So I guess the idea is that step 1 is to look at the Hamiltonian and see if particle number is manifestly conserved or not? (and I suppose if it is conserved, one needs a more complicated method to figure out what the H matrix is?))

• Regarding your addendum, it's pretty easy to enforce particle conservation by only including terms like $c^{\dagger} c$ with equal numbers of creation and annihilation operators. I don't see why this makes writing down the hamiltonian more complicated. Also note that we sometimes consider particle non-conserving hamiltonians for fermions, for example in studying superconductors.
– d_b
Feb 16, 2019 at 1:39

Yes. Each fermion site is a two-dimensional Hilbert space (since the site can only have either zero or one fermion). Denote $$|1\rangle = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$$ as the state with the fermion occupied, and $$|0 \rangle = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$$ as the fermion unoccupied. The creation operator has the action $$c^\dagger |0\rangle = |1\rangle$$, $$c^\dagger |1\rangle = 0$$, so $$c^\dagger = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$$. Likewise, $$c = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}$$. These operators can almost be tensor product-ed together the same way as you do in the spin case. The "almost" is because fermion operators on different sites anti-commute. To enforce that in the second quantization, you need to attach an $$e^{\nu_j}$$, where $$\nu_j = \sum_{i=1}^{j-1} c_i^\dagger c_i = \sum_{i=1}^{j-1} n_i$$, onto the operator at site $$j$$ (see this, for example, on second quantization). This means that $$c_3^\dagger$$ should be represented $$e^{n} \otimes e^n \otimes c^\dagger \otimes I \otimes I \otimes \ldots$$.
Response to addendum: This procedure works fine if the number of particles is conserved. Simply generate the $$2N$$ possible operators (the $$N$$ creation operators and $$N$$ annihilation operators), and then obtain your Hamiltonian via matrix multiplication and addition. However, it brings up an important point, which is relevant if you care about efficiently diagonalizing these sorts of matrices. If our Hamiltonian has some sort of symmetry, then the resulting matrix will be block diagonal. For example, if our Hamiltonian conserves the total number of particles, then any matrix element connecting states with a different number of fermions will be zero. In other words, you can decompose the $$2^N$$ dimensional matrix into a $$1$$ dimensional matrix acting on the zero particle state, an $$N$$ dimensional matrix acting on the $$N$$ possible single particle states, the ($$N$$ choose $$2$$) possible two particle states, and so on, which allows for more efficient diagonalization.