# Bragg-Williams approximation for the Ising model

Again, we consider a system with the following Hamiltonian:

${\displaystyle {\mathcal {H}}=-J\sum _{\left\langle ij\right\rangle }S_{i}S_{j}-\sum _{i}H_{i}S_{i}}$
In this case our random variables ${\displaystyle \varphi }$ are the spins: ${\displaystyle \varphi _{\alpha }=S_{i}}$, and our order parameter is ${\displaystyle \left\langle \varphi _{\alpha }\right\rangle =m_{i}}$ (we will later see that this is exactly the local magnetization also in the mean field approximation). For what we have previously stated we must build a single-particle probability density ${\displaystyle \rho _{i}=f(m_{i})}$ in terms of this order parameter, such that ${\displaystyle \operatorname {Tr} \rho _{i}=1}$ and ${\displaystyle \operatorname {Tr} (\rho _{i}S_{i})=m_{i}}$. For example, since we have two constraints on ${\displaystyle \rho _{i}}$ we could use a linear expression for the probability density involving two parameters[1]:
${\displaystyle \rho _{i}=a\delta _{S_{i},1}+b(1-\delta _{S_{i},1})}$
where ${\displaystyle a}$ and ${\displaystyle b}$ are, respectively, the probability that ${\displaystyle S_{i}=+1}$ and ${\displaystyle S_{i}=-1}$. Inserting this expression of ${\displaystyle \rho _{i}}$ in the two constraints we get:
${\displaystyle \operatorname {Tr} \rho _{i}=a+b=1\quad \qquad \operatorname {Tr} (\rho _{i}S_{i})=\left\langle S_{i}\right\rangle _{\rho _{i}}=a-b=m_{i}}$
and therefore:
${\displaystyle a={\frac {1+m_{i}}{2}}\quad \qquad b={\frac {1-m_{i}}{2}}}$
Thus:
${\displaystyle \rho _{i}={\frac {1+m_{i}}{2}}\delta _{S_{i},1}+{\frac {1-m_{i}}{2}}(1-\delta _{S_{i},1})}$
We are now able to compute the two terms that contribute to the free energy ${\textstyle F_{\rho _{\text{m.f.}}}=\left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}+k_{B}T\sum _{\alpha }\left\langle \ln \rho _{\alpha }\right\rangle _{\rho _{\text{m.f.}}}}$ of the system.

First of all:

${\displaystyle \left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}=-J\sum _{\left\langle ij\right\rangle }\left\langle S_{i}S_{j}\right\rangle _{\rho _{\text{m.f.}}}-\sum _{i}H_{i}\left\langle S_{i}\right\rangle _{\rho _{\text{m.f.}}}}$
However, since ${\textstyle \rho _{\text{m.f.}}=\prod _{i}\rho _{i}}$ (namely the degrees of freedom are independent, as we have already seen), then ${\displaystyle \left\langle S_{i}S_{j}\right\rangle _{\rho _{\text{m.f.}}}=\left\langle S_{i}\right\rangle _{\rho _{i}}\left\langle S_{j}\right\rangle _{\rho _{j}}}$. Furthermore, if ${\displaystyle g(S_{i})}$ is a generic function of ${\displaystyle S_{i}}$ then:
{\displaystyle {\begin{aligned}\left\langle g(S_{i})\right\rangle _{\rho _{\text{m.f.}}}=\operatorname {Tr} (g(S_{i})\rho _{i})=\sum _{S_{i}=\pm 1}g(S_{i})\rho _{i}=\\{}\\=\sum _{S_{i}=\pm 1}g(S_{i})\left({\frac {1+m_{i}}{2}}\delta _{S_{i},1}+{\frac {1-m_{i}}{2}}(1-\delta _{S_{i},1})\right)={\frac {1+m_{i}}{2}}g(1)+{\frac {1-m_{i}}{2}}g(-1)\end{aligned}}}
and so if ${\displaystyle g(S_{i})=S_{i}}$ we have ${\displaystyle \left\langle S_{i}\right\rangle _{\rho _{\text{m.f.}}}=m_{i}}$: our order parameter is (as we expected) the local magnetization of the system. Therefore:
${\displaystyle \left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}=-J\sum _{\left\langle ij\right\rangle }m_{i}m_{j}-\sum _{i}H_{i}m_{i}}$

Then, for the other term of ${\displaystyle F}$ we have:

${\displaystyle \sum _{i}\left\langle \ln \rho _{i}\right\rangle _{\rho _{\text{m.f.}}}=\sum _{i}\operatorname {Tr} (\rho _{i}\ln \rho _{i})=\sum _{i}\left({\frac {1+m_{i}}{2}}\ln {\frac {1+m_{i}}{2}}+{\frac {1-m_{i}}{2}}\ln {\frac {1-m_{i}}{2}}\right)}$

The total energy of the system will be:

{\displaystyle {\begin{aligned}F_{\rho _{\text{m.f.}}}=\left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}+k_{B}T\sum _{i}\left\langle \ln \rho _{i}\right\rangle _{\rho _{\text{m.f.}}}=\\=-J\sum _{\left\langle ij\right\rangle }m_{i}m_{j}-\sum _{i}H_{i}m_{i}+k_{B}T\sum _{i}\left({\frac {1+m_{i}}{2}}\ln {\frac {1+m_{i}}{2}}+{\frac {1-m_{i}}{2}}\ln {\frac {1-m_{i}}{2}}\right)\end{aligned}}}
We now have to minimize ${\displaystyle F_{\rho _{\text{m.f.}}}}$ with respect to ${\displaystyle m_{i}}$:
{\displaystyle {\begin{aligned}0={\frac {\partial F_{\rho _{\text{m.f.}}}}{\partial m_{i}}}=-2J\sum _{j\in {\text{n.n.}}(i)}m_{j}-H_{i}+k_{B}T\left({\frac {1}{2}}\ln {\frac {1+m_{i}}{2}}+{\frac {1}{2}}-{\frac {1}{2}}\ln {\frac {1-m_{i}}{2}}-{\frac {1}{2}}\right)=\\=-2J\sum _{j\in {\text{n.n.}}(i)}m_{j}-H_{i}+{\frac {k_{B}T}{2}}\ln {\frac {1+m_{i}}{1-m_{i}}}\end{aligned}}}
where ${\displaystyle j\in {\text{n.n.}}(i)}$ means that the ${\displaystyle j}$-th site is a nearest neighbour of the ${\displaystyle i}$-th one. Now, recalling that:
${\displaystyle \tanh ^{-1}x={\frac {1}{2}}\ln {\frac {1+x}{1-x}}\qquad |x|<1}$
(and surely ${\displaystyle |m_{i}|<1}$) we can write:
${\displaystyle k_{B}T\tanh ^{-1}m_{i}=2J\sum _{j\in {\text{n.n.}}(i)}m_{j}+H_{i}}$
and inverting the hyperbolic tangent:
${\displaystyle m_{i}=\tanh \left[\beta \left(2J\sum _{j\in {\text{n.n.}}(i)}m_{j}+H_{i}\right)\right]}$
We have again found the self-consistency equation for the magnetization that we have already encountered in the Weiss mean field theory for the Ising model! This is again a confirmation that all mean field theories are equivalent.

We can thus see that in the Bragg-Williams approximation also the ${\displaystyle \beta }$ and ${\displaystyle \alpha }$ exponents are the same of the Weiss mean field theory; in fact, we have seen in Critical exponents of Weiss mean field theory for the Ising model that they come from the expansion of the free energy density for small values of the magnetization when ${\displaystyle H=0}$. If we now set ${\displaystyle H=0}$ into ${\displaystyle F_{\rho _{\text{m.f.}}}}$, so that ${\displaystyle m_{i}=m\quad \forall i}$, we get:

${\displaystyle {\frac {F_{\rho _{\text{m.f.}}}}{N}}=-Jzm^{2}+k_{B}T\left({\frac {1+m}{2}}\ln {\frac {1+m}{2}}+{\frac {1-m}{2}}\ln {\frac {1-m}{2}}\right)}$
If we now expand the logarithm for small ${\displaystyle m}$, we get:
{\displaystyle {\begin{aligned}f(m)=\lim _{N\to \infty }{\frac {F_{\rho _{\text{m.f.}}}}{N}}=-Jzm^{2}+k_{B}T\left(-\ln 2+{\frac {m^{2}}{2}}+{\frac {m^{4}}{12}}+\cdots \right)\sim \\\sim -k_{B}T\ln 2+\left({\frac {k_{B}T}{2}}-Jz\right)m^{2}+{\frac {k_{B}T}{2}}m^{4}\end{aligned}}}
Therefore we see that the behaviour of ${\displaystyle f}$ near ${\displaystyle m=0}$ is the same of equation Weiss mean field theory for the Ising model , and so both the exponents ${\displaystyle \beta }$ and ${\displaystyle \alpha }$ turn out to be equal to what we have determined[2] in Critical exponents of Weiss mean field theory for the Ising model.

We can also say something[3] in the case ${\displaystyle H\neq 0}$. Supposing that our system is uniform, i.e. ${\displaystyle m_{i}=m\quad \forall i}$, we can rewrite ${\displaystyle F_{\rho _{\text{m.f.}}}}$ as:

${\displaystyle {\tilde {f}}:={\frac {F_{\rho _{\text{m.f.}}}}{N}}=-Hm+f_{0}}$
where:
${\displaystyle f_{0}=-Jzm^{2}+k_{B}T\left({\frac {1+m}{2}}\ln {\frac {1+m}{2}}+{\frac {1-m}{2}}\ln {\frac {1-m}{2}}\right)}$
and since we are looking for the absolute minimum of ${\displaystyle {\tilde {f}}}$:
${\displaystyle 0={\frac {\partial {\tilde {f}}}{\partial m}}=-H-f'_{0}=-H-2Jzm+{\frac {k_{B}T}{2}}\ln {\frac {1+m}{1-m}}}$
which gives the self-consistency equation we already know:
${\displaystyle h=-2\beta Jzm+\tanh ^{-1}m}$
As we can see from the figures below the number of the possible solutions depends on the temperature of the system:

Solutions of the self-consistency equation for ${\displaystyle T>T_{c}}$
Solutions of the self-consistency equations for ${\displaystyle T

In particular if ${\displaystyle 2\beta Jz<1}$ (i.e. ${\displaystyle T>T_{c}}$), expanding the right hand side of the self-consistency equation we get a positive linear term (${\displaystyle -2\beta Jzm>0}$), so that the behaviour of ${\displaystyle -\beta Jzm+\tanh ^{-1}m}$ is as shown in the first figure, and the equation has only one solution. On the other hand if ${\displaystyle T then the linear term changes sign and ${\displaystyle -\beta Jzm+\tanh ^{-1}m}$ behaves as in the second figure: in this case if ${\displaystyle |h|}$ is small enough there are three possible solutions, which we have called ${\displaystyle m_{1}}$, ${\displaystyle m_{2}}$ and ${\displaystyle m_{3}}$. These are all extrema of ${\displaystyle {\tilde {f}}}$, but how can we understand which is a minimum or a maximum? And above all, which of them is the global minimum. If we suppose ${\displaystyle h>0}$ to be large there will be only one solution, ${\displaystyle m_{3}}$, and as ${\displaystyle h}$ decreases also ${\displaystyle m_{1}}$ and ${\displaystyle m_{2}}$ will appear; we can therefore argue that for the continuity of ${\displaystyle {\tilde {f}}}$ the solution ${\displaystyle m_{3}}$ is still a minimum also when ${\displaystyle m_{1}}$ and ${\displaystyle m_{2}}$ are present. Similarly, if we take ${\displaystyle h<0}$ we can conclude that also ${\displaystyle m_{1}}$ is a minimum; therefore ${\displaystyle m_{2}}$ will necessarily be a maximum. Now, in order to see which between ${\displaystyle m_{1}}$ and ${\displaystyle m_{3}}$ is the global minimum of ${\displaystyle {\tilde {f}}}$ let us take ${\displaystyle h>0}$ and compute:

${\displaystyle {\tilde {f}}(m_{3})-{\tilde {f}}(m_{1})=\int _{m_{1}}^{m_{3}}{\tilde {f}}'(m)dm=\int _{m_{1}}^{m_{3}}(f'_{0}-h)dm}$
From the second figure we see that this is equal to the area enclosed by the the straight line ${\displaystyle h={\text{const.}}}$ and the graph of ${\displaystyle f'_{0}}$, which is clearly negative if ${\displaystyle h>0}$. Therefore ${\displaystyle {\tilde {f}}(m_{3})<{\tilde {f}}(m_{1})}$ and so we conclude that ${\displaystyle m_{3}}$ is always the global minimum of ${\displaystyle {\tilde {f}}}$. Similarly, when ${\displaystyle h<0}$ the global minimum of ${\displaystyle {\tilde {f}}}$ is ${\displaystyle m_{1}}$.

This means that as soon as ${\displaystyle h}$ changes sign the global minimum of ${\displaystyle {\tilde {f}}}$ changes abruptly from ${\displaystyle m_{3}}$ to ${\displaystyle m_{1}}$. We are thus obtaining the phenomenology that is indeed observed for a magnet when we change the external field ${\displaystyle H}$. In other words the sets of points ${\displaystyle (h(m),m)}$ are exactly the graphs of the ${\displaystyle (H,M)}$ phase diagram we have seen in Phase transitions and phase diagrams, i.e. the graphs of the magnetization seen as a function of the external field.

1. This form of ${\displaystyle \rho _{i}}$ is very general, and does not depend on the fact that the degrees of freedom of the system can only assume two values: if there is a different number of possible states, say ${\displaystyle n}$, then ${\displaystyle \rho _{i}}$ can be written in the same form, but ${\displaystyle a}$ will be the probability of one state while ${\displaystyle b}$ the probability of the remaining ${\displaystyle n-1}$ ones. We will shortly see this when we will apply the Bragg-Williams approximation to the Potts model.
2. Note that also the temperature of the transition is still the same, considering also the ${\displaystyle 1/2}$ factor we have already mentioned.
3. These considerations apply in general also in the other mean field theories considered, but we show them now.