# Hubbard-Stratonovich mean field theory for the Ising model

We have already encountered (see The role of interaction range) the Hubbard-Stratonovich identity with the saddle point approximation to compute the partition function of a one-dimensional Ising model with long-range interactions. The same technique can also be used to formulate mean field theories for systems with short-ranged interactions, and it is one of the most useful ones.

We therefore start from the most general Hamiltonian of a two-dimensional Ising model with nearest-neighbour interactions:

${\displaystyle {\mathcal {H}}=-{\frac {1}{2}}\sum _{\left\langle i,j\right\rangle }J_{ij}S_{i}S_{j}-\sum _{i}H_{i}S_{i}}$
whose partition function is:
${\displaystyle Z_{N}=\operatorname {Tr} e^{\beta \left({\frac {1}{2}}\sum _{\left\langle i,j\right\rangle }J_{ij}S_{i}S_{j}+\sum _{i}H_{i}S_{i}\right)}}$
The idea of Hubbard-Stratonovich mean field theory is the same as in The role of interaction range , namely substituting the quadratic term ${\displaystyle S_{i}S_{j}}$ in the Hamiltonian with a linear one, introducing some auxiliary field over which we integrate. In order to do so we must generalize the Hubbard-Stratonovich identity.

Theorem

Let ${\displaystyle {\boldsymbol {A}}}$ be a real symmetric matrix and ${\displaystyle {\vec {b}}}$ an arbitrary vector. Then:

${\displaystyle {\frac {1}{\sqrt {\det {\boldsymbol {A}}}}}e^{{\frac {1}{2}}b_{i}({\boldsymbol {A}}^{-1})_{ij}b_{j}}=\int _{-\infty }^{+\infty }\prod _{i=1}^{N}{\frac {e^{-{\frac {1}{2}}x_{i}{\boldsymbol {A}}_{ij}x_{j}+x_{i}b_{i}}}{\sqrt {2\pi }}}dx_{i}}$
where for the sake of simplicity we have used Einstein's notation on repeated indices.

We omit the proof of this result.

Therefore, defining:

${\displaystyle b_{i}:=S_{i}\qquad \qquad \;({\boldsymbol {A}}^{-1})_{ij}:=J_{ij}\qquad \qquad \;x_{i}:=\varphi _{i}}$
we get:
${\displaystyle e^{{\frac {1}{2}}J_{ij}S_{i}S_{j}}={\frac {1}{\sqrt {\det J}}}\int _{-\infty }^{+\infty }\prod _{i=1}^{N}{\frac {e^{-{\frac {1}{2}}\varphi _{i}(J^{-1})_{ij}\varphi _{j}+\varphi _{i}S_{i}}}{\sqrt {2\pi }}}d\varphi _{i}}$
and the partition function becomes:
${\displaystyle Z_{N}={\frac {1}{\sqrt {(2\pi )^{N}\det J}}}\int _{-\infty }^{+\infty }\prod _{i=1}^{N}e^{-{\frac {\beta }{2}}\varphi _{i}(J^{-1})_{ij}\varphi _{j}}\operatorname {Tr} _{\lbrace S_{i}\rbrace }e^{\beta (\varphi _{i}+H_{i})S_{i}}d\varphi _{i}}$
where we have explicitly written that the trace is performed over the spin variables ${\displaystyle S_{i}}$. Now that the terms in the exponential are linear in ${\displaystyle S_{i}}$ the partition function is easier to compute; let us change variable defining:
${\displaystyle \psi _{i}=\varphi _{i}+H_{i}}$
so that:
${\displaystyle Z_{N}={\frac {1}{\sqrt {(2\pi )^{N}\det J}}}\int _{-\infty }^{+\infty }\prod _{i=1}^{N}e^{-{\frac {\beta }{2}}(\psi _{i}-H_{i})(J^{-1})_{ij}(\psi _{j}-H_{j})}\operatorname {Tr} _{\lbrace S_{i}\rbrace }e^{\beta \psi _{i}S_{i}}d\psi _{i}}$
However:
${\displaystyle \operatorname {Tr} _{\lbrace S_{i}\rbrace }e^{\beta \psi _{i}S_{i}}=\sum _{\lbrace S_{i}=\pm 1\rbrace }e^{\beta \psi _{i}S_{i}}=\prod _{i}\sum _{S=\pm 1}e^{\beta \psi _{i}S}=\prod _{i}2\cosh(\beta \psi _{i})}$
and since ${\displaystyle x=e^{\ln x}}$ we have:
${\displaystyle \operatorname {Tr} _{\lbrace S_{i}\rbrace }e^{\beta \psi _{i}S_{i}}=e^{\sum _{i}\ln \left[2\cosh(\beta \psi _{i})\right]}}$
Therefore:
${\displaystyle Z_{N}={\frac {1}{\sqrt {(2\pi )^{N}\det J^{-1}}}}\int _{-\infty }^{+\infty }\prod _{i=1}^{N}e^{-\beta {\mathcal {L}}(\psi _{i},H_{i},J_{ij})}d\psi _{i}}$
where:
${\displaystyle {\mathcal {L}}(\psi _{i},H_{i},J_{ij})={\frac {1}{2}}\sum _{i,j}(\psi _{i}-H_{i})(J^{-1})_{ij}(\psi _{j}-H_{j})-{\frac {1}{\beta }}\sum _{i}\ln \left[2\cosh(\beta \psi _{i})\right]}$
Let us note that there is no trace left of the original degrees of freedom ${\displaystyle S_{i}}$ in the partition function, and that it has the form of a functional integral, since in the continuum limit (${\displaystyle N\to \infty }$ and the distance between lattice sites tending to zero) the auxiliary fields become functions of the position, ${\displaystyle \psi ({\vec {r}})}$. Let us also note that until now we still have not done any approximation (we have just rewritten the partition function); this comes into play right now: since the partition function has the form of a functional integral we must find some approximate ways to compute it. The simplest of all possible approximations is the saddle point approximation (see the appendix The saddle point approximation), which essentially consists in approximating the integral with the maximum value of the integrand. In other words we approximate:
{\displaystyle {\begin{aligned}\int _{-\infty }^{+\infty }\prod _{i=1}^{N}e^{-\beta {\mathcal {L}}(\psi _{i},H_{i},J_{ij})}d\psi _{i}\approx e^{-\beta {\mathcal {L}}({\overline {\psi }}_{i},H_{i},J_{ij})}\quad \Rightarrow \\\Rightarrow \quad Z_{N}\approx {\frac {1}{\sqrt {(2\pi )^{N}\det J}}}e^{-\beta {\mathcal {L}}({\overline {\psi }}_{i},H_{i},J_{ij})}\end{aligned}}}
where the stationary solutions are those which satisfy:
${\displaystyle {\frac {\delta {\mathcal {L}}}{\delta \psi _{k}}}_{|{\overline {\psi }}_{k}}=0\qquad \forall k}$
In this case we will have:
${\displaystyle \sum _{j}(J^{-1})_{kj}({\overline {\psi }}_{j}-H_{j})-\tanh(\beta {\overline {\psi }}_{k})=0\qquad \forall k}$
multiplying both sides by ${\displaystyle J_{ik}}$ and summing over ${\displaystyle k}$, taking advantage of the fact that ${\displaystyle J_{ij}=J_{ji}}$ and in the end renaming ${\displaystyle k}$ with ${\displaystyle j}$ we get:
${\displaystyle {\overline {\psi }}_{i}=H_{i}+\sum _{j}J_{ij}\tanh(\beta {\overline {\psi }}_{j})\qquad \forall i}$
Let us see the relation between the auxiliary fields ${\displaystyle {\overline {\psi }}_{i}}$ and the order parameter of the system ${\displaystyle m_{i}}$. We know that:
${\displaystyle m_{k}=-{\frac {\partial F}{\partial H_{k}}}}$
and from the saddle point approximation we have:
${\displaystyle e^{-\beta F}=Z\approx e^{-\beta {\mathcal {L}}({\overline {\psi }}_{i},H_{i},J_{ij})}\quad \Rightarrow \quad F={\mathcal {L}}({\overline {\psi }}_{i},H_{i},J_{ij})}$
Therefore:
{\displaystyle {\begin{aligned}m_{k}=-{\frac {\partial F}{\partial H_{k}}}=-{\frac {\partial }{\partial H_{k}}}{\mathcal {L}}({\overline {\psi }}_{i},H_{i},J_{ij})=\\=-{\frac {\partial }{\partial H_{k}}}\left[{\frac {1}{2}}\sum _{i,j}({\overline {\psi }}_{i}-H_{i})(J^{-1})_{ij}({\overline {\psi }}_{j}-H_{j})-{\frac {1}{\beta }}\sum _{i}\ln \left(2\cosh(\beta {\overline {\psi }}_{i})\right)\right]=\\=\sum _{j}(J^{-1})_{kj}({\overline {\psi }}_{j}-H_{j})=\tanh(\beta {\overline {\psi }}_{k})\end{aligned}}}
Thus:
${\displaystyle {\overline {\psi }}_{i}=k_{B}T\tanh ^{-1}m_{i}\quad \quad \forall i}$
and plugging this into ${\displaystyle {\overline {\psi }}_{i}}$:
${\displaystyle k_{B}T\tanh ^{-1}m_{i}=H_{i}+\sum _{i,j}J_{ij}m_{j}\quad \Rightarrow \quad m_{i}=\tanh \left[\beta \left(H_{i}+\sum _{i,j}J_{ij}m_{j}\right)\right]}$
We can note that this is a more general form of the self-consistency equation for the magnetization that we have found in Weiss mean field theory. In fact, if we set ${\displaystyle H_{i}=H={\text{const.}}}$ (so that also ${\displaystyle m_{i}=m={\text{const.}}}$) and choose a nearest-neighbour interaction (i.e. ${\displaystyle J_{ij}}$ is equal to a constant ${\displaystyle J}$ if the ${\displaystyle i}$-th and ${\displaystyle j}$-th sites are nearest neighbours, otherwise is null), calling ${\displaystyle z}$ the coordination number of the lattice we get exactly what we found in Weiss mean field theory for the Ising model. This means that in this approximation we deduce the same conclusions we have seen there that come from this self-consistency equation; in particular we will have that the temperature of the phase transition when ${\displaystyle H=0}$ is again ${\displaystyle T_{c}=zJ/k_{B}}$ and that the critical exponents ${\displaystyle \delta }$ and ${\displaystyle \gamma }$ are, respectively, 3 and 1. We can thus already see that all the mean field theories are equivalent when applied to the same system (which is something reasonable since the approximation we make is always the same). In the following sections we will see many other mean field theories, applied also to fluids (see Mean field theories for fluids), and we will always find the same values for the critical exponents. This is another expression of the fact that different systems, like an Ising model or a classical fluid, belong to the same universality class (even more, as we have shown in Ising model and fluids they are actually equivalent) and thus behave similarly near a critical point.