Self-consistent mean field algorithm

Mean-field approximation

The full hamiltonian is:

$$\hat{H}=\widehat{H_0}+\hat{V}=\sum _{ij}{h}_{ij}{c}_i^{†}{c}_j+\frac{1}{2}\sum _{ijkl}{v}_{ijkl}{c}_i^{†}{c}_j^{†}{c}_l{c}_k$$

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^{†}=\sum _k{U}_{ik}{b}_k^{†},$$

such that the unperturbed groundstate wavefunction is:

$$|0⟩={\mathrm{\Pi }}_{E_i\le \mu }{b}_i^{†}|\text{vac}⟩.$$

Based on this definition, we define the normal ordering operator

$:ABC...:$ such that it fulfills:

$$:ABC...:|0⟩=0$$

which practically means it orders

$b_i$ operators based on whether its above or below the Fermi level

$\mu$.

Under this definition of normal ordering, we define the Wick's expansion of the interaction term:

$$\begin{array}{c}{c}_i^{†}{c}_j^{†}{c}_l{c}_k=:c_i^{†}{c}_j^{†}{c}_l{c}_k:\hfill \\ +\overline{c_i^{†}{c}_j^{†}}:c_l{c}_k:+\overline{c_i^{†}{c}_k}:c_j^{†}{c}_l:-\overline{c_i^{†}{c}_l}:c_j^{†}{c}_k:+\overline{c_l{c}_k}:c_i^{†}{c}_j^{†}:-\overline{c_j^{†}{c}_k}:c_i^{†}{c}_l:+\overline{c_j^{†}{c}_l}:c_i^{†}{c}_k:\\ \hfill \overline{c_i^{†}{c}_j^{†}}\overline{c_l{c}_k}-\overline{c_i^{†}{c}_l}\overline{c_j^{†}{c}_k}+\overline{c_i^{†}{c}_k}\overline{c_j^{†}{c}_l}\end{array}$$

where the overline defines Wick's contraction:

$$\overline{AB}=AB-:AB:.$$

The expectation value of the interaction with respect to the

$|0⟩$ is:

$$⟨0|c_i^{†}{c}_j^{†}{c}_l{c}_k|0⟩=⟨0|\overline{c_i^{†}{c}_j^{†}}\overline{c_l{c}_k}-\overline{c_i^{†}{c}_l}\overline{c_j^{†}{c}_k}+\overline{c_i^{†}{c}_k}\overline{c_j^{†}{c}_l}|0⟩$$

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:

$$AB\approx ⟨A⟩B+A⟨B⟩-⟨A⟩⟨B⟩$$

onto the contractions such that we get:

$$⟨c_i^{†}{c}_j^{†}{c}_l{c}_k⟩\approx ⟨c_i^{†}{c}_j^{†}⟩⟨c_l{c}_k⟩+⟨c_i^{†}{c}_k⟩⟨c_j^{†}{c}_l⟩-⟨c_i^{†}{c}_l⟩⟨c_j^{†}{c}_k⟩$$

note

$⟨AB⟩\approx ⟨A⟩⟨B⟩$ assuming mean-field.

To consider excitations from the groundstate, we make use of the mean-field approximation defined above:

$$\begin{array}{c}{c}_i^{†}{c}_j^{†}{c}_l{c}_k\approx \hfill \\ ⟨c_i^{†}{c}_j^{†}⟩c_l{c}_k+⟨c_i^{†}{c}_k⟩c_j^{†}{c}_l-⟨c_i^{†}{c}_l⟩c_j^{†}{c}_k+⟨c_l{c}_k⟩c_i^{†}{c}_j^{†}-⟨c_j^{†}{c}_k⟩c_i^{†}{c}_l+⟨c_j^{†}{c}_l⟩c_i^{†}{c}_k+\\ \hfill ⟨c_i^{†}{c}_j^{†}⟩⟨c_l{c}_k⟩+⟨c_i^{†}{c}_k⟩⟨c_j^{†}{c}_l⟩-⟨c_i^{†}{c}_l⟩⟨c_j^{†}{c}_k⟩\end{array}$$

Where we made use of the following operations:

$$:c_i^{†}{c}_j^{†}{c}_l{c}_k:\approx 0$$

$$\overline{c_i^{†}{c}_k}\overline{c_j^{†}{c}_l}\approx ⟨\overline{c_i^{†}{c}_k}⟩\overline{c_j^{†}{c}_l}+\overline{c_i^{†}{c}_k}⟨\overline{c_j^{†}{c}_l}⟩-⟨\overline{c_i^{†}{c}_k}⟩⟨\overline{c_j^{†}{c}_i}⟩=⟨c_i^{†}{c}_k⟩\overline{c_j^{†}{c}_l}+\overline{c_i^{†}{c}_k}⟨c_j^{†}{c}_l⟩-⟨c_i^{†}{c}_k⟩⟨c_j^{†}{c}_l⟩$$

$$\overline{c_i^{†}{c}_k}:c_j^{†}{c}_l:\approx ⟨\overline{c_i^{†}{c}_k}⟩:c_j^{†}{c}_l:+\overline{c_i^{†}{c}_k}⟨:c_j^{†}{c}_l:⟩-⟨\overline{c_i^{†}{c}_k}⟩⟨:c_j^{†}{c}_l:⟩=⟨\overline{c_i^{†}{c}_k}⟩:c_j^{†}{c}_l:$$

$$⟨\overline{c_i^{†}{c}_k}⟩=⟨c_i^{†}{c}_k-:c_i^{†}{c}_k:⟩=⟨c_i^{†}{c}_k⟩$$

Without any superconducting terms, the form simplifies to:

$$\begin{array}{c}{c}_i^{†}{c}_j^{†}{c}_l{c}_k\approx ⟨c_i^{†}{c}_k⟩c_j^{†}{c}_l-⟨c_i^{†}{c}_l⟩c_j^{†}{c}_k-⟨c_j^{†}{c}_k⟩c_i^{†}{c}_l+⟨c_j^{†}{c}_l⟩c_i^{†}{c}_k+\hfill \\ \hfill ⟨c_i^{†}{c}_k⟩⟨c_j^{†}{c}_l⟩-⟨c_i^{†}{c}_l⟩⟨c_j^{†}{c}_k⟩\end{array}$$

Finite size

Coulomb interaction

We simplify the interaction term through the MF approximation to get:

$$V=\frac{1}{2}\sum _{ijkl}{v}_{ijkl}{c}_i^{†}{c}_j^{†}{c}_l{c}_k\approx \frac{1}{2}\sum _{ijkl}{v}_{ijkl}[⟨c_i^{†}{c}_k⟩c_j^{†}{c}_l-⟨c_j^{†}{c}_k⟩c_i^{†}{c}_l-⟨c_i^{†}{c}_l⟩c_j^{†}{c}_k+⟨c_j^{†}{c}_l⟩c_i^{†}{c}_k]$$
(assuming no superconductivity)

and an additional constant part:

$$V_0=\frac{1}{2}\sum _{ijkl}{v}_{ijkl}(⟨c_j^{†}{c}_l⟩⟨c_i^{†}{c}_k⟩-⟨c_j^{†}{c}_k⟩⟨c_i^{†}{c}_l⟩).$$

The interaction reads:

$$\begin{gathered} v_{ijkl}=\iint w_i^{\ast }(r)w_j^{\ast }(r^{\prime })V(r,r^{\prime })w_k(r)w_l(r^{\prime })drd{r}^{\prime }= \\ \iint V(|r-r^{\prime }|)w_i^{\ast }(r)w_k(r)w_j^{\ast }(r^{\prime })w_l(r^{\prime })drd{r}^{\prime } \end{gathered}$$

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:

$$v_{ijkl}^{{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_l}=\iint V(|r-r^{\prime }|)w_{i\times {\sigma }_i}^{\ast }(r)w_{k\times {\sigma }_k}(r)w_{j\times {\sigma }_j}^{\ast }(r^{\prime })w_{l\times {\sigma }_l}(r^{\prime })drd{r}^{\prime }{\delta }_{{\sigma }_i{\sigma }_k}{\delta }_{{\sigma }_j{\sigma }_l}$$

On a fine tight-binding model, we have:

$$v_{ijkl}^{{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_l}=v_{ij}{\delta }_{ik}{\delta }_{jl}{\delta }_{{\sigma }_i{\sigma }_k}{\delta }_{{\sigma }_j{\sigma }_l}$$

where

$v_{ij}=V(r_i,r_j)$.

We shall re-define

$i$ index to absorb spin:

$${\delta }_{ik}\times {\delta }_{{\sigma }_i{\sigma }_k}\to {\delta }_{ik}$$

in this notation the above reads:

$$v_{ijkl}=v_{ij}{\delta }_{ik}{\delta }_{jl}$$

The mean-field terms are:

$$⟨c_i^{†}{c}_j⟩=⟨{\mathrm{\Psi }}_F|c_i^{†}{c}_j|{\mathrm{\Psi }}_F⟩$$

whereas

$|{\mathrm{\Psi }}_F⟩={\mathrm{\Pi }}_{i=0}^{N_F}{b}_i^{†}|0⟩$. 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^{†}=\sum _k{U}_{ik}{b}_k^{†}$$

where

$$U_{ik}=⟨i|{\psi }_k⟩.$$

That gives us:

$$c_i^{†}{c}_j=\sum _{k,l}{U}_{ik}{U}_{lj}^{\ast }{b}_k^{†}{b}_l$$

and its expectation value gives us the mean-field … field

$F_{ij}$:

$$F_{ij}=⟨c_i^{†}{c}_j⟩=\sum _{k,l}{U}_{ik}{U}_{lj}^{\ast }⟨{\mathrm{\Psi }}_F|b_k^{†}{b}_l|{\mathrm{\Psi }}_F⟩=\sum _k{U}_{ik}{U}_{kj}^{\ast }$$

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$:

$$\begin{array}{c}{V}_{nm}=\frac{1}{2}\sum _{ijkl}{v}_{ijkl}⟨n|[⟨c_i^{†}{c}_k⟩c_j^{†}{c}_l-⟨c_j^{†}{c}_k⟩c_i^{†}{c}_l-⟨c_i^{†}{c}_l⟩c_j^{†}{c}_k+⟨c_j^{†}{c}_l⟩c_i^{†}{c}_k]|m⟩=\hfill \\ \frac{1}{2}\sum _{ijkl}{v}_{ijkl}[+{\delta }_{jn}{\delta }_{lm}{F}_{ik}-{\delta }_{in}{\delta }_{lm}{F}_{jk}-{\delta }_{jn}{\delta }_{km}{F}_{il}+{\delta }_{in}{\delta }_{km}{F}_{jl}]=\\ \frac{1}{2}[\sum _{ik}{v}_{inkm}{F}_{ik}-\sum _{jk}{v}_{njkm}{F}_{jk}-\sum _{il}{v}_{inml}{F}_{il}+\sum _{jl}{v}_{njml}{F}_{jl}]=\\ \hfill -\sum _{ij}{F}_{ij}(v_{inmj}-v_{injm})\end{array}$$

where I used the

$v_{ijkl}=v_{jilk}$ symmetry from Coulomb.

For a tight-binding grid, the mean-field correction simplifies to:

$$\begin{array}{c}{V}_{nm}=-\sum _{ij}{F}_{ij}(v_{inmj}-v_{injm})=\hfill \\ -\sum _{ij}{F}_{ij}{v}_{in}{\delta }_{im}{\delta }_{nj}+\sum _{ij}{F}_{ij}{v}_{in}{\delta }_{ij}{\delta }_{nm}=\\ \hfill -F_{mn}{v}_{mn}+\sum _i{F}_{ii}{v}_{in}{\delta }_{nm}\end{array}$$

the first term is the exchange interaction whereas the second one is the direct interaction.

Similarly, the constant offset term reads:

$$\begin{array}{c}{V}_0=\frac{1}{2}\sum _{ijkl}{v}_{ijkl}(F_{jl}{F}_{ik}-F_{jk}{F}_{il})=\hfill \\ \frac{1}{2}\sum _{ijkl}{v}_{ij}{\delta }_{ik}{\delta }_{jl}(F_{jl}{F}_{ik}-F_{jk}{F}_{il})\\ \hfill =\frac{1}{2}\sum _{ij}{v}_{ij}(F_{ii}{F}_{jj}-F_{ji}{F}_{ij})\end{array}$$

where we identify the first term as the exchange (mixes indices) and the right one as the direct (diagonal in indicies).

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_{nm}^{ij}$$

and we will Fourier transform the upper indices into k-space:

$$O_{mn}(k)=\sum _{ij}{O}_{nm}^{ij}{e}^{-ik(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:

$$V_{nm}^{ij}=-F_{mn}^{ij}{v}_{mn}^{ij}+\sum _{r,p}{F}_{pp}^{rr}{v}_{pn}^{rj}{\delta }_{nm}{\delta }^{ij}$$

Lets first consider the second (direct) term. Lets express the corresponding

$F$ term in k-space:

$$F_{pp}^{rr}=\int e^{ik(R_r-R_r)}{F}_{pp}(k)dk=\int F_{pp}(k)dk$$

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:

$$\sum _r{v}_{pn}^{rj}=\int v_{pn}(k)e^{ik{R}_j}\sum _r{e}^{-ik{R}_r}dk=\int v_{pn}(k)e^{ik{R}_j}{\delta }_{k,0}dk=v_{pn}(0)$$

We are finally ready to Fourier transform the main result. Invoking convolution theorem and the results above gives us:

$$V_{nm}(k)=\sum _p{F}_{pp}^0{v}_{pn}(0){\delta }_{nm}-F_{mn}(k)⊛v_{mn}(k)=V_n^D-F_{mn}(k)⊛v_{mn}(k)$$

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.

The constant offset is:

$$\begin{gathered} V_0=\frac{1}{2}\sum _{r,s}{\rho }_r{v}_{rs}(0){\rho }_s- \\ \frac{1}{2}tr[{\int }_{BZ}(F⊛v)(k)F(k)dk] \end{gathered}$$

Short-Range interactions

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:

$$V_{nm}(\mathbf{R})=V_n^D\delta (\mathbf{R})-F_{mn}(\mathbf{R})v_{mn}(\mathbf{R})$$

(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.