Mean field theory for the Blume-Emery-Griffiths model

The Blume-Emery-Griffiths model (often shortened in "BEG" model) is a lattice gas model used to describe the properties of the superfluid transition of ${\displaystyle {}^{4}}$He, which when it is cooled under approximately ${\displaystyle T_{c}=2.7{\text{ K}}}$ undergoes a continuous phase transition[1] from fluid to superfluid; superfluids exhibit several interesting properties, like having zero viscosity. The BEG model is used to describe what happens when we add some ${\displaystyle {}^{3}}$He to the system; it does not consider quantum effects, but only the "messing up" due to the ${\displaystyle {}^{3}}$He impurities. Experimentally when ${\displaystyle {}^{3}}$He is added to ${\displaystyle {}^{4}}$He the temperature of the fluid-superfluid transition decreases. For small concentration of ${\displaystyle {}^{3}}$He the mixture remains homogeneous, and the only effect is the change of ${\displaystyle T_{c}}$; however, when the concentration ${\displaystyle x}$ of ${\displaystyle {}^{3}}$He reaches the critical value ${\displaystyle {\overline {x}}\approx 0.67}$, ${\displaystyle {}^{3}}$He and ${\displaystyle {}^{4}}$He separate into two phases (just like oil separates from water) and the ${\displaystyle \lambda }$ transition becomes first-order (namely, discontinuous). The transition point ${\displaystyle (x_{t},T_{t})}$ where the system shifts from a continuous ${\displaystyle \lambda }$ transition to a discontinuous one is that where the phase separation starts and is called tricritical point. The BEG model was introduced to describe such a situation.

As we have anticipated, it is a lattice gas model and so it is based on an Ising-like Hamiltonian (see Ising model and fluids). On the sites of this lattice we define a variable ${\displaystyle S_{i}}$ which can assume the values ${\displaystyle -1}$, ${\displaystyle 0}$ and ${\displaystyle +1}$: we decide that when an ${\displaystyle {}^{4}}$He atom is present in a lattice site then ${\displaystyle S_{i}=\pm 1}$, while when ${\displaystyle S_{i}=0}$ it means that the site is occupied by an ${\displaystyle {}^{3}}$He atom. We then define our order parameter to be ${\displaystyle m=\left\langle S_{i}\right\rangle }$; in the Ising model ${\displaystyle \left\langle S_{i}^{2}\right\rangle }$ can only be equal to 1, while in this case it can be either 0 or 1: we can thus interpret ${\displaystyle \left\langle S_{i}^{2}\right\rangle }$ as the concentration of ${\displaystyle {}^{4}}$He atoms, and ${\textstyle x:=1-\left\langle S_{i}^{2}\right\rangle }$ as the fraction of ${\displaystyle {}^{3}}$He. We also define ${\displaystyle \Delta :=\mu _{{}^{3}{\text{He}}}-\mu _{{}^{4}{\text{He}}}}$ to be the difference of the chemical potentials of ${\displaystyle {}^{3}}$He and ${\displaystyle {}^{4}}$He; since this parameter is related to the number of ${\displaystyle {}^{3}}$He and ${\displaystyle {}^{4}}$He atoms, we expect that when ${\displaystyle x\to 0}$ (namely, there is only ${\displaystyle {}^{4}}$He) ${\displaystyle \Delta \to -\infty }$, while if ${\displaystyle x\to 1}$ then ${\displaystyle \Delta \to +\infty }$.

We consider the following Hamiltonian for the system[2]:

${\displaystyle {\mathcal {H}}=-J\sum _{\left\langle ij\right\rangle }S_{i}S_{j}+\Delta \sum _{i}S_{i}^{2}-\Delta N}$
where ${\displaystyle N}$ is the total number of lattice sites. Since we want to apply the second variational method that we have seen, we write the mean field probability density as:
${\displaystyle \rho _{\text{m.f.}}=\prod _{i}\rho _{i}=\prod _{i}\rho (S_{i})}$
and the free energy:
${\displaystyle F_{\rho _{\text{m.f.}}}=\left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}+k_{B}T\sum _{i}\operatorname {Tr} (\rho _{i}\ln \rho _{i})}$
The mean value of the Hamiltonian is:
${\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.}}}+\Delta \sum _{i}\left\langle S_{i}^{2}\right\rangle _{\rho _{\text{m.f.}}}-N\Delta }$
and since ${\displaystyle \left\langle S_{i}S_{j}\right\rangle _{\rho _{\text{m.f.}}}=\left\langle S_{i}\right\rangle _{\rho _{\text{m.f.}}}\left\langle S_{j}\right\rangle _{\rho _{\text{m.f.}}}}$ (it's the fundamental hypothesis of mean field theories) we get[3]:
${\displaystyle \left\langle {\mathcal {H}}\right\rangle _{\rho _{\text{m.f.}}}=-{\frac {1}{2}}JNz\left\langle S_{i}\right\rangle _{\rho _{\text{m.f.}}}^{2}+N\Delta \left\langle S_{i}^{2}\right\rangle _{\rho _{\text{m.f.}}}-N\Delta }$
where ${\displaystyle z}$ is the coordination number of the lattice. Therefore, the free energy of the system is:
${\displaystyle F_{\rho _{\text{m.f.}}}=-{\frac {1}{2}}JNz\left(\operatorname {Tr} (\rho _{i}S_{i})\right)^{2}+N\Delta \operatorname {Tr} (\rho _{i}S_{i}^{2})-N\Delta +k_{B}TN\operatorname {Tr} (\rho _{i}\ln \rho _{i})}$
We now must minimize this expression with respect to ${\displaystyle \rho _{i}}$, with the constraint ${\displaystyle \operatorname {Tr} \rho _{i}=1}$. Since we have:

{\displaystyle {\begin{aligned}{\frac {\delta }{\delta \rho _{j}}}\left(\operatorname {Tr} (\rho _{i}S_{i})\right)^{2}=2\operatorname {Tr} (\rho _{i}S_{i})\cdot S_{j}=2\left\langle S_{i}\right\rangle S_{j}=2mS_{j}\\{\frac {\delta }{\delta \rho _{j}}}\operatorname {Tr} (\rho _{i}S_{i}^{2})=S_{j}^{2}\quad \qquad {\frac {\delta }{\delta \rho _{j}}}\operatorname {Tr} (\rho _{i}\ln \rho _{i})=\ln \rho _{j}+1\end{aligned}}}

then:

${\displaystyle {\frac {\delta }{\delta \rho _{i}}}F_{\rho _{\text{m.f.}}}=-JNzmS_{i}+N\Delta S_{i}^{2}+Nk_{B}T\ln \rho _{i}+Nk_{B}T=0}$

${\displaystyle \ln \rho _{i}=\beta JzmS_{i}-\beta \Delta S_{i}^{2}-1\quad \Rightarrow \quad \rho (S_{i})={\frac {1}{A}}e^{\beta (zJmS_{i}-\Delta S_{i}^{2})}}$
where we have reabsorbed ${\displaystyle e^{-1}}$ into the normalization constant ${\displaystyle A}$. From the constraint ${\displaystyle \operatorname {Tr} \rho _{i}=1}$ we find:
${\displaystyle A=1+2e^{-\beta \Delta }\cosh(\beta zJm)}$
Substituting this expression of ${\displaystyle \rho _{i}}$ into ${\displaystyle F_{\rho _{\text{m.f.}}}}$, after some mathematical rearrangement we get:
${\displaystyle {\frac {F_{\rho _{\text{m.f.}}}}{N}}={\frac {z}{2}}Jm^{2}-k_{B}T\ln \left[1+2e^{-\beta \Delta }\cosh(\beta zJm)\right]-\Delta }$
In order to find the equilibrium state for any ${\displaystyle T}$ and ${\displaystyle \Delta }$, we must minimize this expression with respect to ${\displaystyle m}$. If we expand ${\displaystyle F_{\rho _{\text{m.f.}}}}$ for small ${\displaystyle m}$, keeping in mind the Taylor expansions ${\displaystyle \cosh x\sim x+x^{2}/2+O(x^{4})}$ and ${\displaystyle \ln(1+x)\sim 1+x-x^{2}/2+O(x^{3})}$ we get:
${\displaystyle {\frac {F_{\rho _{\text{m.f.}}}}{N}}{\stackrel {m\approx 0}{\sim }}a+{\frac {b}{2}}m^{2}+{\frac {c}{4}}m^{4}+{\frac {d}{6}}m^{6}+\cdots }$
where:
${\displaystyle a=-k_{B}T\ln(1+2e^{-\beta \Delta })-\Delta \qquad b=zJ\left(1-{\frac {zJ}{\delta k_{B}T}}\right)\qquad c={\frac {zJ}{2\delta ^{2}}}(\beta zJ)^{3}\left(1-{\frac {\delta }{3}}\right)}$
and ${\displaystyle \delta =1+e^{\beta \Delta }/2}$; ${\displaystyle d}$ turns out to be always positive. Note that unlike the Ising model in the Weiss approximation (see Weiss mean field theory for the Ising model) in this case both the quadratic and the quartic terms, ${\displaystyle b}$ and ${\displaystyle c}$, can change sign when the parameters assume particular values. Let us also note that the order parameter of the system, namely the concentration of ${\displaystyle {}^{3}}$He, is:
${\displaystyle x=1-\left\langle S_{i}^{2}\right\rangle ={\frac {1}{1+2e^{-\beta \Delta }\cosh(\beta zJm)}}}$

Therefore, in the disordered phase (both ${\displaystyle {}^{3}}$He and ${\displaystyle {}^{4}}$He are present) we have ${\displaystyle m=0}$ and the concentration of ${\displaystyle {}^{3}}$He becomes:

${\displaystyle x={\frac {1}{1+2e^{-\beta \Delta }}}=1-{\frac {1}{\delta }}}$

This way we can determine how the temperature of the ${\displaystyle \lambda }$ transition depends on ${\displaystyle x}$; in fact, the critical temperature will be the one that makes ${\displaystyle b}$ change sign, so we can determine it from the condition ${\displaystyle b=0}$:

${\displaystyle T_{c}={\frac {zJ}{k_{B}\delta }}}$
Since as we have just seen ${\displaystyle 1/\delta =1-x}$, we have:
${\displaystyle T_{c}(x)=T_{c}(0)(1-x)}$
where ${\displaystyle T_{c}(0)=zJ/k_{B}}$. The other transition (from the continuous ${\displaystyle \lambda }$ to the discontinuous one) will occur when the quartic term ${\displaystyle c}$ changes sign, and so we can determine the critical value of ${\displaystyle x}$ at which it occurs from the condition ${\displaystyle c=0}$:
${\displaystyle 1-{\frac {\delta _{t}}{3}}=0\quad \Rightarrow \quad \delta _{t}=3\quad \Rightarrow \quad x_{t}=1-{\frac {1}{3}}={\frac {2}{3}}=0.{\overline {6}}}$
which is in astonishingly good agreement with the experimental result of ${\displaystyle x_{t}\approx 0.67}$.

1. This is generally called ${\displaystyle \lambda }$ transition, because the plot of the specific heat as a function of the temperature has a shape that resembles a ${\displaystyle \lambda }$.
2. We do not justify why it has this form; furthermore, the Hamiltonian we are considering is in reality a simplification of the original one.
3. We are introducing the ${\displaystyle 1/2}$ factor for later convenience. As we have already stated, this is perfectly possible if we change the convention of the sum on nearest neighbours.