Selected Problems | Random physics

$\newcommand{\I}{\mathrm{i}} \newcommand{\E}{\mathrm{e}} \newcommand{\D}{\mathop{}\!\mathrm{d}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\braket}[1]{\langle{#1}\rangle} \newcommand{\bbraket}[1]{\langle\!\langle{#1}\rangle\!\rangle} \newcommand{\bm}{\boldsymbol}$

Lectures on statistical mechanics.

Selected problems: solutions

4 Identical particles

Identity of particles. Find the joint probability density of 2 bosons (and 2 fermions) to be at positions $\bm x_1$ and $\bm x_2$,

$$ \begin{equation} \label{rh} \braket{\bm x_1 \bm x_2 | \rho(H) | \bm x_1 \bm x_2} \end{equation} $$

using the canonical quantum ensemble

$$ \begin{equation} \label{h} \rho(H) = \frac{\E^{-H/T}}{Z_2(T,V)}, \quad H = \frac{\bm p_1^2 + \bm p_2^2}{2m} \end{equation} $$

density matrix; $H$ is the hamiltonian of the two identical particles.

Show that it is given by the formula

$$ \begin{equation} \label{xrx} \braket{\bm x_1 \bm x_2 | \rho(H) | \bm x_1 \bm x_2} = \frac{1}{2\lambda^6 Z_2} \big(1 \pm \E^{-2\pi |\bm x_1 - \bm x_2|^2/\lambda^2}\big) \end{equation} $$

where $\lambda$ is the thermal length at temperature $T$. Deduce the expression of the two particles partition function $Z_2(T,V)$, and compare the result with the classical one.

Solution

The quantum state of the two particles is best described in the momentum representation $\ket{\bm k}$,

$$\bm p \ket{k} = \bm k \ket{k} \quad (\hbar = 1)$$

in which the momentum operator $\bm p$ is diagonal. For bosons we have

$$\ket{k_1k_2}_B = \frac{\ket{k_1k_2} + \ket{k_2k_1}}{\sqrt{2}} = \ket{k_1k_2}_+$$

and

$$\ket{k_1k_2}_F = \frac{\ket{k_1k_2} - \ket{k_2k_1}}{\sqrt{2}} = \ket{k_1k_2}_-$$

for fermions.

In this representation the two particles hamiltonian is diagonal:

$$ \begin{equation} \label{hkk} H = \left( \frac{k_1^2}{2} + \frac{k_2^2}{2} \right) \ket{k_1k_2}\bra{k_1k_2}_\pm \quad (m = 1), \end{equation} $$

and is symmetric in $(k_1, k_2)$. Note that in the present units, the thermal length is $\lambda=(2\pi/T)^{1/2}$ at temperature $T$ ($k_\mathrm{B} = 1$).

Using ($\ref{hkk}$) in ($\ref{rh}$), and introducing the unity decomposition

$$ 1 = \frac{1}{2}\sum_{k_1, k_2} \ket{k_1k_2}\bra{k_1k_2}_\pm $$

where the factor $1/2$ avoids counting twice the same basis projector (their are symmetric in the indices $k_1k_2$), we obtain

$$ \begin{equation} \braket{x_1 x_2 | \E^{-H/T} | x_1x_2} = \frac{1}{2}\sum_{k_1, k_2} \braket{x_1x_2|k_1k_2}_\pm \braket{k_1k_2|x_1x_2}_\pm \exp\left( -\frac{k_1^2}{2T} - \frac{k_2^2}{2T} \right). \end{equation} $$

In the limit of a large volume the summations go to integrals:

$$ \begin{equation} \label{xrxkk} \braket{x_1 x_2 | \E^{-H/T} | x_1x_2} = \frac{V^2}{2} \int \frac{\D \bm k_1}{(2\pi)^3} \frac{\D \bm k_2}{(2\pi)^3} \braket{x_1x_2|k_1k_2}_\pm \braket{k_1k_2|x_1x_2}_\pm \exp\left( -\frac{k_1^2}{2T} - \frac{k_2^2}{2T} \right). \end{equation} $$

and the $\braket{x|k}$ wavefunction becomes:

$$ \braket{x|k} = \frac{\E^{\I \bm k \cdot \bm x}}{\sqrt{V}} $$

which, after substitution into ($\ref{xrxkk}$), gives

$$ \begin{equation} \braket{x_1 x_2 | \E^{-H/T} | x_1x_2} = \frac{1}{4} \int \frac{\D \bm k_1}{(2\pi)^3} \frac{\D \bm k_2}{(2\pi)^3} \Big( 2 \pm \E^{\I \bm k_1 \cdot (\bm x_1 - \bm x_2)} \E^{\I \bm k_2 \cdot (\bm x_2 - \bm x_1)} \pm \E^{\I \bm k_1 \cdot (\bm x_2 - \bm x_1)} \E^{\I \bm k_2 \cdot (\bm x_1 - \bm x_2)} \big)\exp\left( -\frac{k_1^2}{2T} - \frac{k_2^2}{2T} \right). \end{equation} $$

or, noting that the two integrals factorise:

$$ \begin{equation} \braket{x_1 x_2 | \E^{-H/T} | x_1x_2} = \frac{1}{2\lambda^6} \pm \frac{1}{2}\left[\int \frac{\D \bm k}{(2\pi)^3} \E^{\I \bm k \cdot (\bm x_1 - \bm x_2)} \E^{-\frac{k^2}{2T}} \right]^2. \end{equation} $$

The gaussians are readily integrated to give

$$ \begin{equation} \label{xxrxx} \braket{x_1 x_2 | \E^{-H/T} | x_1x_2} = \frac{1}{2\lambda^6} \left( 1 \pm \E^{T|\bm x_1 - \bm x_2|^2} \right) \end{equation} $$

which would coincide with the desired result ($\ref{xrx}$), once the original units are restored.

Integrating now on $x_1x_2$ we deduce the two particle partition function. The first term gives the ideal gas result $(V/\lambda)^2$, the second one contains the entropy corrections resulting form the identity of particles:

$$ Z_2(T,V) = \frac{1}{2!}\left(\frac{V}{\lambda} \right)^2 \pm \frac{V}{2! 2^{3/2}\lambda^3} $$

or equivalently

$$ Z_2(T,V) = \frac{1}{2!}\left(\frac{V}{\lambda} \right)^2 \left(1 \pm \frac{\lambda^3}{2^{3/2} V} \right). $$

Effective entropic interaction between two fermions (FD) and two bosons (BE). At short distances $r$ fermion repulsion is strong (it diverges in $0$, red curve), while bosons attract with an maximum energy of $T \ln2$ (blue curve).

A physical picture of the entropy contribution can be obtained form the observation that ($\ref{xxrxx}$) defines an effective interaction

$$ w(r) = - T \ln \left( 1 \pm \E^{- Tr^2} \right), \quad r = |\bm x_1 - \bm x_2| $$

(note that $\rho \sim \E^{-w/T}$ can be written as $w \sim -T \ln \rho$). The above figure represents the two curves, for bosons ($+$) and fermions ($-$), as a function of the normalized distance $\sqrt{T}r = \sqrt{2\pi}r/\lambda$ (in units of $\sqrt{m/\hbar^2}$).

5

Bose-Einsein condensate. Demonstrate that $N$ free bosons in two dimension do not form a condensate; show in particular that the critical temperature $T_c$ (signaling the quantum dominated regime) vanishes like $T_c \rightarrow 1/\ln N$ in the thermodynamic limit.

In contrast to the free case, demonstrate that in a harmonic trap, a two-dimensional gas of bosons do forms a condensate. Show that in this case the critical temperature is

$$ T_c = \frac{\sqrt{6}\hbar \omega}{\pi} N^{1/2}, $$

where $\omega$ the harmonic oscillator frequency.

Solution

The interesting quantity is the number of particles, and in particular their distribution between the ground and excited energy levels

$$ \begin{equation} \label{Nsum} N = \sum_s \frac{1}{\E^{(\epsilon_s - \mu)/T}-1} \end{equation} $$

For a free gas confined in a two-dimensional box of area $L^2$, we start computing the density of states, in order to transform ($\ref{Nsum}$) into an integral:

$$ \nu(\epsilon) = \frac{\D}{\D \epsilon} \Big[ \frac{L^2}{(2\pi \hbar)^2} \pi p^2\Big] = \frac{L^2}{2\pi} \frac{m}{\hbar^2} = \frac{L^2}{2\pi} $$

where, in the last equality we put $\hbar=m=1$, and we used that $\epsilon_{\bm p} = p^2/2m$. Note also that the ground state energy

$$ \epsilon_0 = \frac{\pi^2\hbar^2}{2mL^2} $$

tends to zero in the thermodynamic limit ($L \rightarrow \infty$). Just in the thermodynamic limit we write

$$ \begin{equation} \label{Nint} N = \frac{L^2}{2\pi} \int_{0}^\infty \frac{\D \epsilon}{\E^{-\mu/T} \E^{\epsilon/T} - 1} = \frac{L^2 T}{2\pi} \ln\big( 1 - \E^{\mu/T} \big) \end{equation} $$

and may conclude that the chemical potential vinishes with the temperature,

$$ \mu = T \ln\big(1 - \E^{2\pi n/T} \big), \quad n = N/L^2 $$

and therefore, there is no transition to a condensate in two dimensions.

To refine the above argument, which depends on taking the thermodynamic limit, let us discuss the distribution of the particles taking $N$ finite and $\mu = 0$ (the critical value for which the integral in Eq. $\ref{Nint}$ diverges)

$$ N - N_0 = \frac{L^2}{2\pi} \int_{\epsilon_0}^\infty \frac{\D \epsilon}{\E^{\epsilon/T} - 1} \approx \frac{L^2 T}{2\pi} \ln \frac{T_c}{\epsilon_0} $$

This equation defines the “critical” temperature $T_c$ as the one for which $N_0 = 0$:

$$ N = \frac{L^2}{2\pi} \int_{\epsilon_0}^\infty \frac{\D \epsilon}{\E^{\epsilon/T_c} - 1} \approx \frac{L^2 T_c}{2\pi} \ln \frac{T_c}{\epsilon_0}. $$

This is a transcendental equation for $T_c$ that can be solved iteratively (for large $N$):

$$ \begin{equation} T_c \approx \frac{2\pi N}{L^2 \ln N}. \end{equation} $$

where we replaced inside the logarithm $T_c/\epsilon_0 = 2L^2 T_c/\pi^2$ by $4N/\pi$, and neglected irrelevant terms. We finally obtain

$$ \begin{equation} T_c \approx \frac{2\pi \hbar^2 n}{m \ln N}, \end{equation} $$

in the original units. This result shows how the “critical” temperature vanishes in the thermodynamic limit. Because of the slow logarithmic dependence on $N$, in the region $T<T_c(N)$ we may observe a fraction of particles in the ground state.

If instead of free particles we have trapped bosons in an harmonic potential of frequency $\omega$, the one-particle energy is (in two dimensions)

$$ \epsilon_n = \hbar \omega(n_x + n_y + 1). $$

Equation ($\ref{Nsum}$) becomes in this case

$$ \begin{equation} \label{Nsumho} N = \sum_{n_x, n_y} \frac{1}{\exp\big(\frac{\hbar \omega n}{T} - \frac{\mu - \hbar \omega}{T}\big)-1}, \quad n = n_x + n_y \end{equation} $$

which shows that the critical chemical potential is $\mu = \hbar \omega$, and the sum is over $(n_x,n_y)$ such that $n = n_x + n_y = \epsilon/\hbar \omega -1$. The density of states, in terms of $n$, is given by

$$ \nu(n) = \frac{\D}{\D n} \frac{n^2}{2} = n. $$

The critical temperature is determined by the condition $\mu = \hbar \omega$. It is convenient to split the sum ($\ref{Nsumho}$) into one term containing the number of particles in the ground state, and one term containing the number of particles in the excited states,

$$ N(T) = N_0(T) + N_e(T), \quad \mu = \hbar \omega $$

where, for the excited states, we can replace the sum by an integral:

$$ N_e(T) = \int_0^\infty \frac{n \D n}{\E^{\hbar\omega n/T} - 1} = \frac{T^2}{(\hbar\omega)^2} \Gamma(2) \zeta(2) = \frac{\pi^2 T^2}{6(\hbar\omega)^2}, $$

from which we obtain, putting $N_0 = 0$ and then $N_e = N$, the value of the critical temperature,

$$ \begin{equation} \label{Tcho} T_c = \hbar \omega\frac{\sqrt{6N}}{\pi}. \end{equation} $$

Therefore, the number of particles in the ground state, for $T<T_c$ is given by

$$ N - N_0 =\frac{\pi^2 T^2}{6(\hbar\omega)^2}, $$

$$ \frac{N_0}{N} = 1 - \left( \frac{T}{T_c} \right)^2. $$

One may think that the critical temperature ($\ref{Tcho}$) is not well defined in the thermodynamic limit (when $N \rightarrow \infty$). In fact, we must consider that $\omega = \omega(N)$. Indeed, the potential energy of the oscillator can be written in terms of the characteristic length of the trap $R$, which determines the density of the bosons trapped in the area $R^2$, and its characteristic energy $U_0$:

$$ \frac{1}{2} m\omega x^2 = \frac{U_0}{2} \left(\frac{x}{R}\right)^2, \quad \omega = \sqrt{\frac{U_0}{m}}\frac{1}{R}, $$

thus, $\omega \sim 1/R$. Using these parameters, the critical temperature writes

$$ T_c = \hbar \sqrt{\frac{U_0}{m}} \frac{1}{R} \frac{\sqrt{6N}}{\pi} = \hbar \sqrt{\frac{6U_0}{\pi^2 m}} \sqrt{n}, \quad n = \frac{N}{R^2} $$

which is indeed an intensive variable. Equivalently, we have $\omega \sim 1/\sqrt{N}$ with $\omega \sqrt{N} = \mathrm{const}$, in the thermodynamic limit.

6 Thermionic emission

Electrons in a metal at temperature $T$ can escape the surface energy barrier $-W$ if their velocity $p_z/m$ (perpendicular to the surface) satisfies $p_z > \sqrt{2 m W}$ ($m$ is the electron mass and $\bm p$ its momentum). Using the Fermi-Dirac distribution $n_\epsilon$ show that the thermionic current is

$$j = \frac{emT^2}{2\pi^2 \hbar^3} \exp\left(-\frac{W - \epsilon_F}{T}\right)$$

where $\epsilon_F$ is the fermi energy ($k_\mathrm{B} = 1$).

Use units such that $e = W = m = \hbar = 1$. Assume that typical particle energies are larger than $W$ (otherwise emission is impossible). Compare the quantum result with the classical one.

Solution

With $k_\mathrm{B} = e = W = m = \hbar =1$, mass, length and time dimensions are:

$$ \begin{gather} \mathrm{M} = m, \\ \mathrm{L} = \sqrt{\frac{\hbar^2}{mW}}, \\ \mathrm{T} = \frac{\hbar}{W}, \end{gather} $$

and energy and electric current density are:

$$ \begin{gather} \mathrm{E} = W, \\ \mathrm{j} = [env] = e \frac{mW^2}{\hbar^3}. \end{gather} $$

Moreover, temperature as usual is measured in energy units (here $W$).

The electron current $j$ is obtained from the formula

$$ \begin{equation} j = 2 \int_{\sqrt{2}}^\infty \frac{\D p_z}{2\pi} \int_{-\infty}^\infty \frac{\D p_x}{2\pi} \int_{-\infty}^\infty \frac{\D p_y}{2\pi} p_z n_{\bm p} \end{equation} $$

where $n_{\bm p}$ is the fermi distribution

$$ n_{\bm p} = \frac{1}{\E^{(\epsilon _{\bm p} - \mu)/T} + 1},\quad \epsilon_{\bm p} = \frac{p^2}{2}$$

$\bm p$ is the electron momentum, and, in particular $p_z/m$, its perpendicular to the surface velocity. The front factor 2 comes from the two spin states of the electron. Note that only particles with $p_z > (2mW)^{1/2}$ contribute to the current. Our task reduces to the computation of the above integral.

The integral over $p_x, p_z$ is elementary in polar coordinates:

$$ \begin{equation} j = 2 \int_{\sqrt{2}}^\infty \frac{\D p_z}{2\pi} p_z \int_0^\infty \frac{\D p p}{2\pi} \frac{1}{\E^{(p_z^2/2 - \mu)/T} \E^{-p^2/2T} + 1}. \end{equation} $$

Indeed, after changing the integration variable to $x = p^2/2T$ we are left to calculate an integral of the form

$$\int_0^\infty \frac{\D x}{a\E^x +1} = \ln(1 + 1/a).$$

The result is

$$ j = \frac{T^2}{2\pi^2} \int_{\sqrt{2}}^\infty \D p_z p_z \ln\big(1 + \E^{-(p_z^2/2 -\mu)/T}\big). $$

This is not elementary, but the physical relevant regime is the one in which the exponential factor is small with respect to the unity: the energies of interest are those satisfying $\epsilon_z = p_z^2/2 > W$, well above the fermi energy $\mu = \epsilon_F$ (knowing that $T < \epsilon_F$ in most metals at laboratory temperatures). In this case, we expand $\ln(1+x) \approx x$ to get

$$ j = \frac{T^2}{2\pi^2} \int_{\sqrt{2}}^\infty \D p_z p_z \E^{-(p_z^2/2 -\mu)/T} = \frac{T^2}{2\pi^2} \E^{-(1-\epsilon_F)/T}. $$

In the original units (we multiply by $\mathrm{j}$) we finally obtain the thermionic current at temperature $T$:

$$ \begin{equation} j = \frac{emT^2}{2\pi^2\hbar^3} \E^{-(W-\epsilon_F)/T}. \end{equation} $$

The corresponding classical result $j_c$, can be obtained by replacing fermi chemical potential by the ideal gas one:

$$\E^{\mu/T} = \frac{N\lambda^3}{2V}$$

which leads to

$$ j_c = e n_0 \sqrt{\frac{T}{2\pi m}} \E^{-W/T} $$

($n_0 = N/V$) a simple numerical estimation shows that this result is wrong. Another way to estimate classically this current is to use the boltzmann factor $n = n_0 \exp(-W/T)$ and the (one-dimensional) thermal speed $v = \sqrt{T/m}$:

$$ j_c = env \sim en_0 \sqrt{\frac{T}{m}} \E^{-W/T} $$

which coincides with the previous result up to a constant numerical factor.

8

Mean field. The hamiltonian of an ising model in a regular lattice of $N$ sites with $q$ neighbors is

$$ H = -J\sum_{\braket{x,y}} s_x s_y $$

where $s_x = \pm 1$ is the spin on site $x$; $y \in N(x)$ is a neighbor of $x$ belonging to the set $N(x)$ of size $q$; $J$ is the exchange parameter. In the Bethe approximation, which is a generalization of the standard mean-field model, one considers a given spin $s_0$ and its neighbors:

$$ H(s_0) = -Js_0 \sum_{y \in N(0)}s_y - h \sum_{y \in N(0)}s_y = H_c, $$

where the last term contains the information about the local field $h$ created at the site $0$ by the neighboring spins $y \in N(0)$. $H_c$ is the “cluster” hamiltonian, it depends on the unknown parameter $h$, which should be determined self-consistently to match the expectation value of the observables (like the magnetization $m = \braket{s_x}$).

Using the “cluster” partition function

$$ Z_c = \mathrm{tr} \, \E^{-H_c(h)/T} = \sum_{s_x = \pm 1} \E^{-H_c(h)/T}, $$

compute $\braket{s_0}$ and $s_x$. From the condition of homogeneity (at equilibrium) $m = \braket{s_0} = \braket{s_x}$ for all $x$, deduce the self-consistent condition

$$ \frac{\cosh^{q-1}[(J+h)/T]}{\cosh^{q-1}[(J-h)/T]} = \E^{2h/T} $$

which determines $T_c$:

$$ T_c \approx \frac{2J}{\ln\frac{q}{q-2}}. $$

Compare the result for $q=4$ (square lattice) with the Onsager exact value of $T_c$.

Show that near the ferromagnetic transition, the magnetization behaves as $m \sim |T-T_c|^{1/2}$:

$$ m^2 \approx \frac{3 q}{q-1} \frac{q}{2} \ln \left( \frac{q}{q-2} \right) \left( 1 - \frac{T}{T_c} \right). $$

Verify that this formula reduces to the standard one (mean field) in the limit $q \rightarrow \infty$.

Solution

The idea of the bethe approximation is to combine the mean-field method, which takes into account the influence of the mean magnetization on each spin, and the near neighbors interactions with a given spin, Our task is to compute self-consistently the “cluster” field associated to a particular site, using translation invariance.

Bethe cluster of site $s_0$ (black dot), $J$ is the exchange energy; the neighbors $s_y$ (gray dots) are responsible of the cluster field $h$.

The partition function of one cluster $Z_c$ (the one associated to $s_0$, for example) is

$$ Z_c = Z(s_0) = \sum_{s_0 = \pm1} \prod_y \Big(\sum_{s_y=\pm1}\Big) \E^{-H(s_0)/T} = \sum_{s_0 = \pm1} \prod_{y=1}^{q} \sum_{s_y=\pm1} \E^{(J/T)s_0s_y + (h/T)s_y}. $$

The sums over the cluster configurations give

$$ Z_c = \left(2 \cosh \frac{J+h}{T} \right)^{q} + \left(2 \cosh \frac{J-h}{T} \right)^{q}. $$

The computation of mean value of the reference spin is straightforward,

$$ \begin{equation} \label{e:s0} \braket{s_0} = \frac{1}{Z_c} \left[\left(2 \cosh \frac{J+h}{T} \right)^{q} - \left(2 \cosh \frac{J-h}{T} \right)^{q}\right]. \end{equation} $$

For the neighbor $x$ of $s_0$ in the cluster we write

$$ \braket{s_x} = \frac{1}{Z_c} \left[ \prod_{y=1}^{q} \sum_{s_y=\pm1} s_x \E^{(J+h)s_y/T} + \prod_{y=1}^{q} \sum_{s_y=\pm1} s_x \E^{(J-h)s_y/T} \right], $$

separating the case $x=y$ to the others we obtain

$$ \braket{s_x} = \frac{1}{Z_c} \left[ 2\sinh \frac{J+h}{T} \left(2 \cosh \frac{J+h}{T} \right)^{q-1} - 2\sinh \frac{J-h}{T} \left(2 \cosh \frac{J-h}{T} \right)^{q-1} \right] $$

Translation invariance imposes that the mean value of a spin do not depend on the site (the equilibrium state is homogeneous) $\braket{s_0} = \braket{s_x}$, therefore, from the previous results we deduce

$$ \begin{equation} \label{e:hc} \frac{\cosh\frac{J+h}{T}}{\cosh\frac{J-h}{T}} = \E^{2h/(q-1)T}, \end{equation} $$

which is a self-consistent equation that determines the possible values of the cluster field $h$ (see figure below). Simple algebra shows that $h$ is given by the roots of the function

$$ f(h,T) = h - \frac{T(q-1)}{2J} \ln \frac{\cosh\frac{J+h}{T}}{\cosh\frac{J-h}{T}} $$

represented on the figure.

For $T>T_c$ only the root at $h=0$ is present; for $T<T_c$ two new roots appear, given the possible values of the self-consistent cluster field. We take $J=1$ and $T=T_c/1.3$, and obtain $h = \pm 2.22$ (dots).

We write a python code to compute the value of $h$ (the dots in the figure):

import numpy as np
from scipy.optimize import root_scalar
def f(h):
    q = 4
    hmax = 1.6
    bc = 1/2 * np.log(q/(q-2)) # inverse critical temperature
    b = 1.3*bc # below the critical temperature
    return h - (q-1)/(2*b) *np.log(np.cosh(b*(1+h)) / np.cosh(b*(1-h)))
sol = root_scalar(f, bracket=[1, 3], method='brentq')

The critical temperature $T_c$ corresponds to the value of $T$ for which $h=0$ becomes an unstable root of $f$. It can be determined form the condition that the slope at $h=0$ of the left hand side of ($\ref{e:hc}$) must equals its right hand side:

$$ \frac{\partial}{\partial h}\left. \frac{\cosh\frac{J+h}{T}}{\cosh\frac{J-h}{T}} \right|_{h=0} = \frac{\partial}{\partial h}\left. \E^{2h/(q-1)T} \right|_{h=0} , $$

$$ 2\tanh \frac{J}{T_c} = \frac{2}{q-1}, $$

or equivalently

$$ \frac{J}{T} = \frac{1}{2} \ln \frac{q}{q-2}. $$

For the square lattice this gives $T_c \approx 2.885$, that you can compare to the exact value of Onsager $2.269$, a huge improvement with respect to the simple mean-field calculation which gives $T_c = 4$ (in units of $J$). To get the previous expression you may show that

$$ 2x = \ln \frac{1+t}{1-t}, \quad t = \tanh x. $$

Now we compute the magnetization $m$ near the critical temperature, when $T \approx T_c$ and $h \approx h_c$. Equation ($\ref{e:s0}$) gives us the relation between the cluster field $h$ and $m = \braket{s_0}$

$$ m = \tanh \frac{qh}{(q-1)T} \approx \frac{qh}{(q-1)T} , $$

where, in the last equality, we used the fact that $h/T \ll 1$ (close to the transition). From the right hand side of the same equation ($\ref{e:s0}$), we obtain

$$ \frac{h}{T} \approx (q-1) \tanh(J/T) \frac{h}{T} - \frac{q-1}{3} \frac{\tanh(J/T)}{\cosh^2(J/T)} \left(\frac{h}{T}\right)^3 + \cdots $$

where we used an expansion in a taylor series. Below $T_c$ we have

$$ (q-1)\tanh(J/T) > 1 = (q-1)\tanh(J/T_c), $$

therefore,

$$ \tanh(J/T) - \tanh(J/T_c) \approx \frac{1}{3} \tanh(J/T_c)\mathrm{sech}^2(J/T_c) \left(\frac{h}{T}\right)^3, $$

where, within the same approximation, we replaced $T=T_c$. Expanding the difference of tangents

$$\tanh(x+\epsilon) - \tanh(x) \approx \mathrm{sech}^2(x) \epsilon,$$

we finally obtain

$$ h = J\sqrt{3(q-1)} \left(1-\frac{T}{T_c}\right)^{1/2} $$

or, in terms of the magnetization per site,

$$ m \sim \left(1-\frac{T}{T_c}\right)^{1/2}, $$

the usual power low with exponent $\beta=1/2$ of the mean field theory. The bethe approximation improves the estimation of the critical temperature, but do not change the critical exponents at the phase transition, showing that these exponents are rather insensitive to the influence of the field created by the near neighbors. Only long range correlation can modify the mean-field results.

9

Mean field and Landau theory. Helium ($^4\mathrm{He}$) liquefies at $T_l \approx 4\,\mathrm{K}$ (at normal pressure), and undergoes a phase transition around $T_c \approx 2.2\,\mathrm{K}$ towards a superfluid state, the so-called $\lambda$-transition. If the helium liquid contains a fraction $c_3$ of the $^3\mathrm{He}$ isotope, the nature of the superfluid transition changes when $x > x_t = c_3/(c_3 + c_4) \approx 0.67$, the tricritical point ($c_4$ is the $^4\mathrm{He}$ concentration): it is second order for $x<x_t$, and first order when the mixture when $x > x_t$, which results in the separation into rich in helium-3 and rich in helium-4 phases. A simple model to describe the properties of the tricritical point $(x_t, T_t)$, was proposed by Blume, Emery and Griffiths (1971). It is based on the hamiltonian

$$ H = -J \sum_{\braket{i,j}} S_i S_j + \Delta \sum_i (S_i^2 - 1), \quad S_i = 0, \pm 1 $$

where $S_i = \pm 1$ represents the helium-4 interacting atoms and $S_i = 0$ the helium-3 atoms, treated as impurities, $\Delta$ controls the number of helium-3 atoms (it is proportional to the difference of chemical potentials); the sums extend to the total number of lattice sites $N$, and $\braket{i,j}$ labels neighbors. The concentration of $^3\mathrm{He}$ is then $x = 1 - \braket{S_i^2}$.

Show that the probability distribution P(S) is given by

$$ P(S) = \frac{\E^{(qJmS - \Delta S^2)/T}}{Z}, \quad Z = 1 + 2\E^{-\Delta/T} \cosh\big(qJm/T\big). $$

You may use the minimization of the free energy $F$, assuming that $S_i$ are independent random variables (mean field approximation):

$$ F = \braket{H} + NT\sum_S P(S) \ln P(S), $$

where we also exploit homogeneity. Compute the expected value of $S^2$:

$$ \braket{S^2} = \frac{2\E^{-\Delta/T} \cosh(qJM/T)}{Z}. $$

Deduce the free energy per atom $f(m; T, \Delta) = F/N$

$$ f = \frac{1}{2} qJm^2 - T \ln \left( 1 + 2 \E^{-\Delta/T}\cosh \frac{qJm}{T} \right) - \Delta. $$

Just above the transition (disordered state) $m = 0$; you can then expand $f(m)$ in powers of $m$ up to the 6th order. Analyse the possible phase transitions using Landau approach. Compute, in particular, the values of the concentration $x_t$ and the temperature $T_t$ at the tricritical point.

Numerical application: plot the phase diagram in the $(x,T)$ plane.

Solution

One essential feature of the mean-field approximation is the neglect of correlation: $\braket{S_xS_y} \sim \braket{S_x}\braket{S_y}$. On the level of the probability distribution function $P(S_1,S_2,\ldots,S_N)$ of the configuration $\mathcal{S} = \{S_1S_2\cdots S_N \mid S_x=0,\pm1,\; \forall x\}$, this assumption means that it factorizes as

$$ P(\mathcal{S}) = \prod_{x=1}^N P(S_x) = P(S)^N $$

where in the last equality we used translation invariance $P(S_x) = P(S)$. As a consequence, the computation of the mean energy is straightforward:

$$ \braket{H} = -\frac{qNJ}{2} \braket{S}^2 + \Delta N \braket{S^2} - \Delta N. $$

Indeed, for instance the $J$ term can be written as

$$ \braket{\sum_{\braket{x,y}} S_x S_y} = \frac{1}{2} \sum_x \sum_{y \in N(x)} \braket{S_x S_y} \approx \frac{1}{2} \sum_x \sum_{y \in N(x)} \braket{S} \braket{S}, $$

which leads to the first term in $\braket{H}$, where $q$ is the number of neighbors and $N(x)$ the neighbors of $x$ set.

The free energy $F = \braket{H} - T S_G$ where $S_G$ is the Gibbs entropy, is then

$$ \frac{F}{N} = -\frac{qJ}{2} \Big[\sum_S P(S) S\Big]^2 + \Delta \sum_S P(S)S^2 - \Delta + T \sum_S P(S) \ln P(S). $$

To find $P(S)$ we minimize $F = F[P(s)]$ subject to the constraint $\sum_S P(S) = 1$:

$$ \frac{1}{N} \frac{\delta}{\delta P(S)} F = -qJ \Big[\sum_S P(S) S\Big] S + \Delta S^2 + T \ln P(S) + (T - \mu) = 0, $$

where $\mu$ is the lagrange multiplier associated to the normalization of $P$. We hence obtain,

$$ P(S) = \frac{1}{Z} \exp\Big[\frac{qJ}{T} m S - \frac{\Delta}{T} S^2\Big], \quad m = \sum_S P(S) S. $$

We note that this is a self-consistent equation ($m$ depends on $P$). To ensure its validity, the mean value of $S$ must satisfy

$$ m = \frac{1}{Z} \Big[\E^{qJm/T - \Delta/T} - \E^{-qJm/T - \Delta/T}\Big] $$

or, taking into account that

$$ Z = 1 + \Big[\E^{qJm/T - \Delta/T} + \E^{-qJm/T - \Delta/T}\Big] $$

we write

$$ \begin{equation} \label{begm} m = \frac{2\E^{-\Delta/T} \sinh(qJm/T)}{1 + 2\E^{-\Delta/T} \cosh(qJm/T)}. \end{equation} $$

Similarly, the computation of $\braket{S^2}$, gives

$$ \begin{equation} \label{begx} \braket{S^2} = \frac{2\E^{-\Delta/T} \cosh(qJm/T)}{1 + 2\E^{-\Delta/T} \cosh(qJm/T)} = 1- x. \end{equation} $$

Using these results, we can write the free energy as a function of the macroscopic variable $m$:

$$ \frac{F(m, T, \Delta)}{N} = f(m, T, \Delta) = -\frac{qJ}{2} m^2 + \frac{2\Delta \E^{-\Delta/T}}{Z} \cosh(qJm/T) - \Delta - TS_G, $$

and

$$ S_G = -\frac{1}{Z} \ln \frac{1}{Z} - \frac{1}{Z} \exp\Big[\frac{qJ}{T} m - \frac{\Delta}{T} \Big]\Big(\frac{qJm}{T} - \frac{\Delta}{T} + \ln \frac{1}{Z}) - \frac{1}{Z} \exp\Big[-\frac{qJ}{T} m - \frac{\Delta}{T} \Big]\Big(-\frac{qJm}{T} - \frac{\Delta}{T} + \ln \frac{1}{Z}) $$

$$ S_G = -\frac{1}{Z} \big[1 + \E^{-\Delta/T} \cosh(qJm/T)\big] - \frac{2qJm}{Z} \E^{-\Delta/T} \sinh(qJm/T) + \frac{2\Delta}{TZ} \E^{-\Delta/T} \cosh(qJm/T) , $$

where we recognize the expressions of $\braket{S}=m$ and $\braket{S^2}$. Using these formulas, we can rewrite the free energy in the form,

$$ \begin{equation} \label{BEGfmt} f(m, T, \Delta) = \frac{qJ}{2}m^2 - T \ln\Big[1 + 2\E^{-\Delta/T} \cosh(qJm/T)\Big] - \Delta. \end{equation} $$

We are now in a position to compute the thermodynamics of the Blume, Emery and Griffiths model.

Possible phase transitions can be monitored by $m$ and $x$, the two order parameters of the helium mixture: $m$ describes the normal-superfluid transition (it vanishes in the normal phase), while $x$ describes the influence of the $^3\mathrm{He}$ (it vanishes in the pure $^4\mathrm{He}$ phase). Note that the free energy can be written as

$$ f = \frac{qJ}{2}m^2 + T \ln(x) - \Delta. $$

We are interested in the line separating the $m=0$ and $m\ne 0$ phases as a function of the $^3\mathrm{He}$ concentration. From ($\ref{begm}$) and ($\ref{begx}$) we get

$$ x_0 = \frac{\E^{\Delta/T}}{\E^{\Delta/T} + 2}, \quad \text{or } \Delta = T \ln \frac{2x_0}{1-x_0},\quad \text{for } m=0, $$

and, dividing the two equations,

$$ x(m) = 1 - m \coth(qJm/T), \quad \text{for } m \ne 0. $$

Numerically, we can then obtain $m = m(T)$ from ($\ref{begm}$) in which we replace $x$ by $x(m)$. The whole phase diagram $(\Delta, T)$, giving the equilibrium values of $(m,x)$, can be obtained by solving ($\ref{begm}$) and ($\ref{begx}$) for each pair of values of $\Delta$ and $T$. In experiments the control parameters are $(x,T)$, hence it is convenient to transform the plane $(\Delta,T)$ into $(x,T)$.

In the case $x=0$ we must recover the usual second order phase transition between the normal and superfluid phases of liquid helium. Expanding $f = f_0 + f_2 m^2 + f_4 m^4$ around $m=0$ we obtain (see the python code below where we use sympy to compute the expansion coefficients $f_n,\; n = 0, 2, 4, 6$)

$$ \begin{gather} f_0 = -\Delta + T\ln\frac{\delta - 1}{\delta} \\ f_2 = \frac{(qJ)^2}{2T\delta} \Big(\frac{T\delta}{qJ} - 1\Big) \\ f_4 = \frac{(qJ)^4}{8T^3 \delta^2} \Big(1 - \frac{\delta}{3}\Big) , \end{gather} $$

where

$$ \delta = \frac{1}{2} \E^{\Delta/T} + 1 = \frac{1}{1-x}. $$

We remark that the second order coefficient $f_2$ becomes negative, the condition for a second order transition, at $T=T_c(x) = qJ/\delta$ or

$$ T_c(x) = qJ(1 - x), $$

which is then the critical temperature of the superfluid transition, it linearly decreases with the helium 3 concentration.

import sympy as sy
from sympy import series
m, T, q, J, d, delta = sy.symbols("m, T, q, J, Delta, delta", real=True)
f = q*J*m**2/2 - T*sy.log(1 + 2*sy.exp(-d/T)*sy.cosh(q*J*m/T)) - d
f_poly = series(f, m, 0, 7)
# simplify each coefficient
f_coeff = []
for n in [0,2,4,6]:
    f_coeff.append(f_poly.coeff(m, n=i).subs(sy.exp(d/T),2*delta-2).simplify())

Moreover, the fourth order coefficient $f_4$, which must be positive for the second order phase transition, becomes negative when $\delta > \delta_c = 3$. Therefore

$$ x_0 = 2/3, \quad \Delta_c = \frac{1}{\delta_c} \ln(2\delta_c -2) = 0.462, \quad T_c = \frac{1}{3} $$

should correspond to the transformation of the second order phase transition into a first order one.

Free energy as a function ($f+\Delta$) of the order parameter $m$, for different values of the helium 3 concentration $x$ and temperature. The critical curve $x_c,T_c$ (gray) shows a flat minimum. For $x<x_c$ the phase transition is of second order, its separates the superfluid to the normal fluid phases. For $x > x_c$ we see instead that the shape of $f$ changes, developing secondary minima at a finite distance from zero: $m = m(T)$ becomes a discontinuous function, signaling the presence of a first order phase transition (red).

In order to investigate the neighborhood of the tricritical point $(x_c,T_c)$ in which superfluid, normal, and separated phases coexist, we must allow $m \ne 0$. This can be done by introducing a fictive field $H$ and to define the Gibbs free energy by $g = f + Hm = g(H)$. This corresponds to adding $-mH$ to the free energy $f$. To expand $g$ in powers of $m$ we must find the function $g = g(H(m))$, where $H$ is considered a function of $m$ (equation of state):

$$ H = \frac{\partial g}{\partial m} = 2am + 4bm^3 + 6cm^5 = H(m), $$

where we wrote the expansion of $g$ in powers of $m$

$$ g = g_0(\Delta,T) + a(\Delta,T) m^2 + b(\Delta,T) m^4 + c(\Delta,T) m^6. $$

This can be done by inverting the function $m=m(H)$ using the self-consistent equation

$$ \begin{equation} m = \frac{2\E^{-\Delta/T} \sinh(qJm/T + H/T)}{1 + 2\E^{-\Delta/T} \cosh(qJm/T + H/T)}. \end{equation} $$

which gives $m$ for $H \ne 0$, and computing the coefficients $(a,b,c)$ from its derivatives with respect to $m$ at $m=0$. For example, the derivative of the previous equation with respect to $m$ gives an equation for $\partial H/\partial m$:

$$2a= \left. \frac{\partial H}{\partial m} \right|_{m=0} \dots$$

A sympy code is shown below.

import sympy as sy
from sympy import series
m, T, q, J, d, delta, H = sy.symbols("m, T, q, J, Delta, delta, H", real=True)
a, b, c = sy.symbols("a, b, c", real=True)
# free energy f(m,H)
g = q*J*m**2/2 - T*sy.log(1 + 2*sy.exp(-d/T)*sy.cosh(q*J*m/T + H/T)) - d
# magnetization
g1 = sy.expand(sy.diff(g, m)/J/q)
# expand to 6th order
g1h = g1.subs(H, 2*a*m + 4*b*m**3 + 6*c*m**5)
# simplification
def subs_m0(f):
    return f.subs(m, 0).subs(sy.exp(d/T), 2*delta-2).simplify()
# m^2 coefficient a
g2 = subs_m0(sy.diff(g1h, m))
a_sol = sy.solve(g2, a)[0]
# m^4 coefficient b
g4 = subs_m0(sy.diff(g1h, m, 3))
b_sol = sy.solve(g4, b)[0].subs(a, a_sol).simplify()
# m^6 coefficient c
g6 = subs_m0(sy.diff(g1h, m, 5))
c_sol = sy.solve(g6, c)[0].subs({a: a_sol, b: b_sol}).simplify()

We get

$$ \begin{gather} a = \frac{qJ}{2} \Big(\frac{T\delta}{qJ} - 1\Big)\\ b = \frac{T\delta^2}{8} \Big(1 - \frac{\delta}{3} \Big)\\ c = \frac{T\delta^3}{12}\Big( \frac{3}{20} \delta^2 - \frac{3}{4}\delta + 1\Big) \end{gather} $$

This is the Landau expansion of the gibbs free energy; note that the conditions to find the tricritical point are as before (where we used $f$), but now the 6th order coefficient is always positive (it was not the case in the previous calculation, which failed already at the forth order, just because of the presence of the tricritical point).

In the region where $x>x_c$ (and $T \sim T_c$) the shape of $g$ (or equivalently of f ($\ref{BEGfmt}$)) shows three separated minima (red line in the figure). This corresponds to a first order transition. The separation line between the normal fluid and the separated mixture in the phase space $(x,T)$ can be determined by the condition that the polynomial $g(m)$ should write

$$ g(m) = g_0 + c m^2 (m^2 - m^2_\pm)^2, $$

in order to have the minima just at $m = 0$ and $m = m_\pm$ ($m_- = - m_+$) with the same value of $g=g_0$. This implies the relation

$$ b^2 = 4 ac. $$

between $a,b,c$, which determines a curve in the $(T,x)$ plane. Therefore, on this curve $T = T(x)$ (for x > x_c$), $a>0$ ($c>0$), then $T> qJ/\delta$:

$$ T_1(\delta) = \frac{24}{\delta} \frac{3\delta^2 - 15\delta + 20}{67\delta^2 - 330\delta + 435}. $$

Near the critical point we have $\delta = 1/(1-x)$, and $T(x)$ is slightly above $1-x$ of the second order transition at lower concentration.

Further information can be obtained in the original paper:

Blume M, Emery VJ, and Griffiths RB, Ising Model for the $\lambda$ Transition and Phase Separation in He3-He4 Mixtures. Phys. Rev. A 4, 1071 (1971) https://doi.org/10.1103/PhysRevA.4.1071.