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

${\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:
$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 $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 ${\boldsymbol {A}}$ be a real symmetric matrix and ${\vec {b}}$ an arbitrary vector. Then:

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

$b_{i}:=S_{i}\qquad \qquad \;({\boldsymbol {A}}^{-1})_{ij}:=J_{ij}\qquad \qquad \;x_{i}:=\varphi _{i}$ we get:
$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:
$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 $S_{i}$ . Now that the terms in the exponential are linear in $S_{i}$ the partition function is easier to compute; let us change variable defining:
$\psi _{i}=\varphi _{i}+H_{i}$ so that:
$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:
$\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 $x=e^{\ln x}$ we have:
$\operatorname {Tr} _{\lbrace S_{i}\rbrace }e^{\beta \psi _{i}S_{i}}=e^{\sum _{i}\ln \left[2\cosh(\beta \psi _{i})\right]}$ Therefore:
$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:
${\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 $S_{i}$ in the partition function, and that it has the form of a functional integral, since in the continuum limit ($N\to \infty$ and the distance between lattice sites tending to zero) the auxiliary fields become functions of the position, $\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:
{\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:
${\frac {\delta {\mathcal {L}}}{\delta \psi _{k}}}_{|{\overline {\psi }}_{k}}=0\qquad \forall k$ In this case we will have:
$\sum _{j}(J^{-1})_{kj}({\overline {\psi }}_{j}-H_{j})-\tanh(\beta {\overline {\psi }}_{k})=0\qquad \forall k$ multiplying both sides by $J_{ik}$ and summing over $k$ , taking advantage of the fact that $J_{ij}=J_{ji}$ and in the end renaming $k$ with $j$ we get:
${\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 ${\overline {\psi }}_{i}$ and the order parameter of the system $m_{i}$ . We know that:
$m_{k}=-{\frac {\partial F}{\partial H_{k}}}$ and from the saddle point approximation we have:
$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:
{\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:
${\overline {\psi }}_{i}=k_{B}T\tanh ^{-1}m_{i}\quad \quad \forall i$ and plugging this into ${\overline {\psi }}_{i}$ :
$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 $H_{i}=H={\text{const.}}$ (so that also $m_{i}=m={\text{const.}}$ ) and choose a nearest-neighbour interaction (i.e. $J_{ij}$ is equal to a constant $J$ if the $i$ -th and $j$ -th sites are nearest neighbours, otherwise is null), calling $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 $H=0$ is again $T_{c}=zJ/k_{B}$ and that the critical exponents $\delta$ and $\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.