* If we tried to solve $N$ particles with $M$ distinct states exactly, the full Hilbert space is $M^N$. Mean-field maps into a much smaller non-interacting system of size $M N$
* Quantitative accuracy is either hit or miss, but gives a good understanding of the qualitative behaviour of the system.
where $s_i$ is the spin of the $i$th particle, $J$ is the exchange coupling, and $h$ is the external magnetic field. The sum is over nearest neighbours.
#### Identifying and applying mean field
It is an interacting problem and thus a headache to solve. Rather than consider all the particles, lets consider a single particle in an **effective field** of other spins. To do that, lets define the average spin $m_{i}\equiv \langle s_{i} \rangle$ and re-write the $s_i s_j$ product into:
where for simplicity I also assumed that the average spin $m_i$ is constant for all particles.
#### Self-consistency
The average mean-field $m$ here acts as a variable, or better yet, **an initial guess**. After solving $H_{\text{MF}}$, we need to make sure that the $m$ is **self-consistent** with the new $s_i$:
$$
m = \frac{1}{N} \sum_i^N \langle s_i \rangle.
$$
So we re-calculate $m$, plug it back into the equation of $H_{\text{MF}}$, and repeat until $m$ converges.
#### Summary
1. Identify the mean-field variables and construct the mean-field Hamiltonian.
2. Guess the initial mean-field.
3. Self-consistency loop:
1. Solve the mean-field Hamiltonian $H_{\text{MF}}$ for the given mean-field.
2. Calculate the new mean-field.
3. Check convergence. If not converged, go back to step 3.1. with the new mean-field.
4. ???
5. Profit.
## Why waste time with this?
Interacting systems have become quite hot research field in condensed matter physics (I mean, take a look at graphene). However, numerical packages to solve them on tight-binding systems lack the following:
* Not many well-maintained packages.
* Code is needlessly complex and documentation is lacking.
* Lack generality.
## Idea behind our implementation
### Identifying mean-fields
#### Real Space
You can find the whole theory [here](https://hackmd.io/@-DUiWUyjQXei-EsdckYO-w/HyEbQhIjo).
Here the the main points. A general particle number preserving interaction with all mean-fields can be written as a sort of [Wick's contraction](https://en.wikipedia.org/wiki/Wick%27s_theorem):
$$
V = \frac{1}{2}\sum_{ijkl} v_{ijkl} c_i^{\dagger} c_j^{\dagger} c_l c_k
whereas $|\Psi_F \rangle = \Pi_{i=0}^{N_F} b_i^\dagger |0\rangle$. To make sense of things, we need to transform between $c_i$ basis (position + internal dof basis) into the $b_i$ basis (eigenfunction of a given mean-field Hamiltonian):
$$
c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger
$$
where $U$ is the matrix of eigenvectors of the mean-field Hamiltonian:
where the second term is the Direct Coulomb interaction, and the first term is the exchange interaction.
#### k-space or translational invariance case
The above works for a finite sized tight-binding model, but what about a periodic system? In that case, we can use the Fourier transform to write the mean-field Hamiltonian in k-space. I'll spare you the details, the final result reads:
where $\rho_{p}$ is the particle density at unit cell site $p$, averaged over a k-grid:
$$
\rho_{p} = \int F_{pp}(k) dk
$$
Once again, the first term (exchange) is purely responsible for the hopping whereas the second term (direct) is a potential term coming from the mean-field.
## Scaling of the algorithm
Lets say $M$ is the number of degrees of freedom within the unit cell, and $N$ is the number of k-points along a direction (we consider 2D problem here). Then the following steps limit the scaling of the algorithm:
* Eigenvalue problem for each k-point: $O(N^2 M^3)$.
* Convolution in k-space: $O(N^4 M^2)$. In this case, this is the most expensive step.