where $|\textrm{vac}\rangle$ is the vacuum state, $b_i^\dagger$ is the creation operator of eigenstate $i$ with energy $E_i$ of $\hat{H}_\text{mf}$ and $\mu$ is the Fermi level.
We relate the $b_i^\dagger$ particles to our original basis via a unitary transformation:
$$
c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger.
$$
## Normal ordering and contractions
Before proceeding, we define the *normal ordering* operation, $:ABC...:$, as a sorting of operators such that all the creation operators $b_i^\dagger$ below the Fermi level are to the left of the annihilation $b_i$ operators below the Fermi levels, and vice versa for the operators above the Fermi level.
Whenever the normal ordered operators acts on the ground state, it gives zero:
$$
:ABC...: | 0 \rangle = 0.
$$
Lastly, we define the *contraction* of two operators $A$ and $B$ as:
$$
\overline{AB} = \hat{A}\hat{B} - :AB:.
$$
## Expansion of the interaction term
We utilize Wick's theorem to expand the interaction term:
where the first term in the second equality is zero due to the normal ordering operation and the second term is zero since we assume deviations from the mean-field ground state are small.
Next, we apply the mean-field approximation to the second term in the expansion:
We assume the dominant part of the ground state wavefunction comes from $\hat{H}_0$. Let's assume $b_i$ operators diagonalize the unperturbed hamiltonian:
$$
c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger,
$$
such that the unperturbed groundstate wavefunction is:
where we can forget about all the normal ordered states since those give zero acting on the unperturbed groundstate. To evaluate this further, we utilize the mean-field approximation:
$$
A B \approx \langle A \rangle B + A \langle B \rangle - \langle A \rangle \langle B \rangle
note $\langle A B \rangle \approx \langle A \rangle \langle B \rangle$ assuming mean-field.
where $c_i^\dagger$ and $c_i$ are creation and annihilation operators respectively for fermion in state $i$.
The first term $\hat{H}_0$ is the non-interacting Hamiltonian which by itself is straightforward to solve in a single-particle basis by direct diagonalizations made easy through packages such as `kwant`.
The second term $\hat{V}$ is the interaction term between two particles.
In order to solve the interacting problem exactly, one needs to diagonalize the full Hamiltonian $\hat{H}$ in the many-particle basis which grows exponentially with the number of particles.
Such a task is often infeasible for large systems and one often needs to resort to approximations.
To consider excitations from the groundstate, we make use of the mean-field approximation defined above:
\iint V(|r - r'|) w_{i}^*(r)w_{k}(r) w_{j}^*(r') w_l(r') dr dr'
$$
whereas $w_i$ is a wannier function on site i (and corresponding dof). Whenever one interchanges $i \to j, k \to l$, the Coulomb term is preserved $v_{ijkl} = v_{jilk}$
To make things more understandable, we are also going to explicitly split up position and spin indices: $i \to i \times \sigma$. In this notation, the Coulomb integral reads:
where operator $\delta \hat{A}$ describes the fluctuations around the order parameter.
Let us consider an additional operator $\hat{B}$ and say we are interested in the product of the two operators $\hat{A}\hat{B}$.
If we assume that the fluctuations $\delta$ are small, we can approximate the product of operators into a sum of single operators and the product of the expectation values:
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):
where we make use of Wicks theorem to simplify the expression and neglect the superconducting pairing and constant offset terms.
$$
c_i^\dagger = \sum_{k} U_{ik} b_k^\dagger
$$
where
:::{admonition} Derivation of the mean-field Hamiltonian with Wicks theorem
:class: dropdown info
```{include} mf_details.md
```
:::
The expectation value terms $\langle c_i^\dagger c_j \rangle$ are the density matrix elements of the ground state of the system:
whereas I assumed the indices for wavefunctions $k,l$ are ordered in terms of increasing eigenvalue. We pop that into the definition of the mean-field correction $V$:
where I used the $v_{ijkl} = v_{jilk}$ symmetry from Coulomb.
For density-density interactions (like Coulomb repulsion) the interaction tensor reads:
where we identify the first term as the exchange (mixes indices) and the right one as the direct (diagonal in indices).
## Translational Invariance
The above assumed a finite tight-binding model - all $nm$-indices contain the position of all atoms (among other dof). In this section tho we want to consider an infinite system with translational invariance.
To begin with we deconstruct a general matrix $O_{nm}$ into the cell degrees of freedom ($nm$) and the position of the the cell itself ($ij$):
$$
O_{nm} \to O^{ij}_{nm}
$$
and we will Fourier transform the upper indices into k-space:
$$
O_{mn}(k) = \sum_{ij} O_{nm}^{ij} e^{-i k (R_i-R_j)}
$$
where I assumed $O$ (and thus all operators I will consider here) is local and thus diagonal in k-space.
Now lets rewrite our main result in the previous section using our new notation:
Notice that in the final expression, there is no $rr$ dependence and thus this term is cell-periodic. Therefore, we shall redefine it as cell electron density $\rho$:
$$
F_{pp}^0 = F_{pp}(R = 0) = \int F_{pp}(k) dk
$$
Now since $\rho$ has no $r$ dependence, we can proceed with the sum:
which does make sense. The first term (direct) is a potential term coming from the mean-field and the second term (exchange) is purely responsible for the hopping.
In the case of short-range interactions, it is much more convenient to go back to real space to both store objects and perform the operations. In real space the mean-field part of the Hamiltonian reads:
(the first term might need some prefactor from Fourier transformation)
where $\mathbf{R}$ is a sequence of integers representing real-space unit cell indices.
### Proposed Algorithm
Given an initial Hamiltonian $H_0 (R)$ and the interaction term $v(R)$ in real-space, the mean-field algorithm is the following:
0. Start with a mean-field guess in real-space: $V(R)$.
1. Fourier transform tight-binding model and the mean-field in real space to a given k-grid: $H_0(R) + V(R) \to H_0(k) + V(k)$
2. Diagonalize and evaluate the density matrix: $H_0(k) + V(k) \to F(k)$
3. Fourier transform the density matrix back to real-space: $F(k) \to F(R)$
4. Evaluate the new mean-field Hamiltonian $V(R)$ according to the equation above.
5. Evaluate self-consistency metric (could be based either on $V(R/k)$ or $F(R/k)$). Based on that, either stop or go back to 1 with a modified $V(R)$ starting guess.
