点击跳转到中文版（Chinese ver.)

The microcanonical ensemble of the previous chapter describes an isolated system: the energy is given once and for all, and every average has to be computed inside a thin energy shell. This is conceptually clean but practically painful, because the shell restriction enters every integral and the entropy $S(E,V)$ has to be obtained before anything else can be computed.

In this chapter we relax the constraint on the energy. The system is put in contact with a heat reservoir, so that its energy fluctuates, and the ensemble average turns into an unrestricted integral over the whole phase space. Then we relax the constraint on the particle number as well, which leads to the grand canonical ensemble and to the chemical potential.

Canonical ensemble

Canonical ensemble

We consider an isolated system made up of two subsystems

1: $N_1, V_1, H_1$
2: $N_2, V_2, H_2$

with

and assume that the interaction between the two subsystems is negligibly small.

As discussed in the previous chapter, a “most probable” decomposition $E = \bar{E}_1 + \bar{E}_2$ of the energy dominates. We are interested in the special case $N_2 \gg N_1$ and $\bar{E}_2 \gg \bar{E}_1$, so that subsystem 2 is a large environment of subsystem 1.

Let $E = E_1 + E_2$ be any decomposition of the energy, fluctuating around the most probable one. The probability of finding system 1 at $(p_1,q_1)$ and system 2 at $(p_2,q_2)$ is $\propto \dd p_1 \dd q_1 \dd p_2 \dd q_2$.

If we look at only system 1, we should sum up all possible microscopic states of system 2 with $E_2 = E - E_1$,

That is, the probability of finding system 1 in a microscopic state within the volume $\dd p_1 \dd q_1$ at $(p_1,q_1)$ is $\propto \dd p_1 \dd q_1\, \Gamma_2(E-E_1)$, and therefore the density distribution is

Since $E_1$ fluctuates only around $\bar{E}_1$, and $\bar{E}_1 \ll E$,

Here $T$ is the temperature of the whole system, and therefore also the temperature of system 1, and $S_2(E)$ is irrelevant to system 1.

Question: why not expand $\Gamma_2(E-E_1)$ itself?

Answer: $\Gamma_2$ is the volume of phase space, and it grows exponentially with $N_2$. Expanding such a function in powers of $E_1$ produces coefficients that grow with $N_2$, and the series is useless beyond the first term. The quantity that is additive and slowly varying is not $\Gamma_2$ but its logarithm, i.e. the entropy, and $\log \Gamma_2$ is extensive. Besides, $\Gamma_2$ is not even smooth on the scale of $\Delta$: it is a coarse-grained volume, defined only up to the shell width. The expansion of $k\log\Gamma_2$ is legitimate because the second-order term is smaller than the first by a factor of order $E_1/E$, and we have assumed $E_1 \ll E$.

Omitting the subscript 1, this leads to

This is the probability distribution of the so-called canonical ensemble. Similar to that of the microcanonical ensemble, $\rho(p,q)$ of the canonical ensemble depends on $(p,q)$ also only through the Hamiltonian $H(p,q)$.

The macroscopic conditions for the canonical ensemble:

• $N$ is fixed;
• the system is in contact with a heat reservoir of temperature $T$, and $T$ is fixed.

In other words, the energy is no longer a macroscopic parameter of the system: it is replaced by the temperature of the reservoir. Fig. 3.1(a) shows the arrangement.

1
\begin{tikzpicture}[>=stealth, thick, font=\small]
2
  % ---- (a) canonical ensemble ----
3
  \draw (0,0) rectangle (2,2);
4
  \node at (1,1.35) {System 1};
5
  \node at (1,0.75) {$N_1,V_1$};
6

7
  \draw (2.9,0) rectangle (6.7,2.8);
8
  \node at (4.8,2.2) {Heat reservoir 2};
9
  \node at (4.8,1.6) {$N_2\gg N_1$};
10
  \node at (4.8,1.05) {$\bar{E}_2\gg\bar{E}_1$};
11
  \node at (4.8,0.5) {$T$ fixed};
12

13
  \draw[<->,red] (2.05,1.0) -- (2.85,1.0);
14
  \node[red,above] at (2.45,1.05) {$E$};
15
  \node at (1,-0.5) {(a)};
16

17
  % ---- (b) grand canonical ensemble ----
18
  \draw (0,-5.2) rectangle (6.7,-1.9);
19
  \node[anchor=north west] at (0.15,-2.0) {Reservoir: $T$, $\mu$};
20

21
  \draw (2.5,-4.3) rectangle (4.9,-2.35);
22
  \fill (2.9,-2.7) circle (2pt);
23
  \fill (4.4,-2.75) circle (2pt);
24
  \fill (3.0,-3.9) circle (2pt);
25
  \fill (4.5,-3.95) circle (2pt);
26
  \node at (3.7,-4.6) {open system: $V$ fixed, $N$ varies};
27

28
  \fill (0.6,-2.5) circle (2pt);
29
  \fill (1.1,-3.4) circle (2pt);
30
  \fill (1.7,-4.3) circle (2pt);
31
  \fill (0.6,-4.7) circle (2pt);
32
  \fill (5.9,-2.4) circle (2pt);
33
  \fill (6.2,-3.4) circle (2pt);
34
  \fill (5.5,-4.3) circle (2pt);
35
  \fill (6.2,-4.7) circle (2pt);
36

37
  \draw[->,red] (2.0,-3.25) -- (2.4,-3.25);
38
  \node[red,above left] at (2.3,-3.3) {$\Delta N$};
39
  \draw[->,red] (5.0,-3.25) -- (5.4,-3.25);
40
  \node[red,above right] at (5.1,-3.3) {$\Delta N$};
41

42
  \node at (1,-5.65) {(b)};
43
\end{tikzpicture}

We define the partition function as the effective volume in $\Gamma$ space, i.e. the effective number of the microscopic states,

then the thermodynamics is obtained from

where $A(V,T)$ is the free energy. To show this, one should prove

(a) $A$ is extensive;
(b) $A$ is related to $U$ by $A = U - TS$ with $S \equiv -(\partial A/\partial T)_V$.

Proof of (b): differentiating the following identity with respect to $\beta = 1/kT$,

we obtain, keeping in mind that $H(p,q)$ does not depend on $\beta$,

Since

where we have used $T = 1/(k\beta)$ and therefore $\partial T/\partial \beta = -T/\beta$. Therefore the bracket is just $A - H + TS$, and averaging gives $A - U + TS = 0$, i.e. $A = U - TS$, with $U = \langle H \rangle$ the internal energy.

Proof of (a): the partition function of a system composed of two independent parts factorizes, $Q_{N_1+N_2} = Q_{N_1}Q_{N_2}$, up to the same binomial counting and the same neglect of $O(\log N)$ terms as in the proof of the extensiveness of $S$ in the previous chapter. Hence $-A/kT = \log Q_N$ is additive, and $A$ is extensive.

Please note that now the temperature is already defined from the beginning: it is the temperature of the reservoir, and it is an independent variable. Then the free energy $A(V,T)$ is a characteristic function, and one may straightforwardly derive

The free energy $A(V,T)$ in the canonical ensemble plays the role of the entropy $S(E,V)$ in the microcanonical ensemble.

In principle, the integration of $\dd p \dd q$ cannot extend to the entire $\Gamma$ space, since from the beginning $\rho(p_1,q_1)$ requires $E_1 \le E$. But we show in the next subsection that the contribution from the region $H \ge E$ is negligibly small. In fact, only one value of the energy dominantly contributes to the partition function. Therefore, practically the integration may extend over all energies. This provides much convenience for calculations in the canonical ensemble.

In addition, the free energy $A(V,T)$ in the canonical ensemble can be explicitly computed. In contrast, the internal energy $U(S,V)$ in the microcanonical ensemble could be potentially obtained only from the computation of $S(E,V)$, and such a solution is usually much more difficult.

The definition of the partition function $Q_N(V,T)$ in the canonical ensemble is similar to the total volume $\Gamma(E)$ of the microscopic states in the microcanonical ensemble. In fact

Question: why and how?

Hint: only the microscopic states with the most probable energy $\bar{E}$, i.e. with the average energy $U = \bar{E}$, contribute to $Q_N(V,T)$. Therefore $Q_N(V,T) = e^{-\bar{E}/kT}\Gamma(\bar{E})$, and taking the logarithm gives $k\log Q_N = -U/T + S$ directly.

A simple argument for the canonical ensemble: suppose that a large isolated system consists of many identical quasi-independent subsystems, and each subsystem is described by the Hamiltonian $H(p,q)$. According to the result of the previous chapter, i.e. the Maxwell-Boltzmann distribution, the probability distribution of finding a subsystem at $(p,q)$ is

The same result therefore follows without repeating the reservoir argument, provided that the interaction between the subsystems is negligible.

Energy fluctuation

Questions:

• In what sense are the microcanonical and canonical ensembles equivalent?
• Why do we need both ensembles?

Any macroscopic observable can be calculated by

and for the energy,

which we may write as

Differentiating both sides with respect to $\beta$, and noting that $\partial U/\partial\beta$ is the only term in which the derivative acts on something other than the exponential, we obtain

Further, since

the bracket is $A - H + TS = U - H$, and therefore $\partial U/\partial\beta + \langle (U-H)^2 \rangle = 0$, i.e.

or, in terms of the deviation $\Delta E = H - \langle H \rangle$,

In fact, this is a special case of the well-known fluctuation-dissipation theorem, which relates the fluctuation of the energy to the response to the external change. Since the left-hand side is a variance, it also follows that $C_V \ge 0$ in equilibrium.

Now, because $\langle H \rangle \sim N$ and $C_V \sim N$, we have, as $N \to \infty$,

as sketched in Fig. 3.2. Therefore, the microcanonical and canonical ensembles are equivalent, in the sense that the energy only fluctuates around a fixed value in both cases. In the canonical ensemble, however, the energy is not given and needs to be computed. Further, this theorem also shows that the contribution of the states with a very large energy is negligible: the weight $e^{-\beta H}$ falls off exponentially, and the fluctuation is small.

1
\begin{tikzpicture}[>=stealth, thick, font=\small, x=0.8cm, y=2.8cm]
2
  \draw[->] (0,0) -- (10.5,0) node[below] {$N$};
3
  \draw[->] (0,0) -- (0,1.05) node[above] {$\frac{\sqrt{\langle(\Delta E)^2\rangle}}{U}$};
4
  \node[below left] at (0,0) {$O$};
5

6
  \draw[blue] plot[domain=0.95:10, samples=120] (\x, {0.85/sqrt(\x)});
7

8
  \draw[dashed] (1,0) -- (1,0.85);
9
  \draw[dashed] (4,0) -- (4,0.425);
10
  \draw[dashed] (4,0.425) -- (0,0.425);
11
  \draw[dashed] (1,0.85) -- (0,0.85);
12

13
  \fill (1,0.85) circle (2pt);
14
  \fill (4,0.425) circle (2pt);
15
  \fill (9,0.2833) circle (2pt);
16

17
  \node[right, xshift=4pt] at (1,0.85) {$N$};
18
  \node[right, xshift=4pt] at (4,0.425) {$4N$};
19
  \node[above] at (9,0.32) {$O(N^{-1/2})$};
20
\end{tikzpicture}

In the figure the relative fluctuation of the energy is drawn against $N$; the curve is $O(N^{-1/2})$, so quadrupling the system halves the relative fluctuation.

Boltzmann statistical theory

Consider $N$ identical quasi-independent particles,

where now the index $i$ labels particles. Then

where $Z(\beta,V) = \sum_{p,q} e^{-\beta\epsilon(p,q)}$ is the partition function of a single particle, the same $Z(\beta,y)$ introduced in the previous chapter, with $y$ replaced by $V$.

Using $\log N! \approx N\log N - N$, we then have $\log Q_N \approx N\log Z - N\log N + N = N\log(Z/N) + N$, i.e.

The constant term $1$ comes from Stirling’s approximation, and it cannot be dropped when the entropy is computed. Then

since $\partial\beta/\partial T = -\beta/T$.

Comparing $Z(\beta,V)$ here with $Z(\beta,y)$ introduced in the previous chapter, the only difference is that the energy for calculating $Z(\beta,V)$ here is not bounded by the total energy. But this difference is negligibly small, since $\epsilon \ll E \propto N$.

Some applications of canonical ensemble

Equipartition theorem

Let the Hamiltonian be

then $\left\langle \frac{1}{2} a_i p_i^2 \right\rangle = \frac{1}{2}kT$.

Proof:

because the integrations over all the other variables cancel between the numerator and the denominator. For the numerator, integrating by parts,

where the boundary term vanishes because $H \to \infty$ as $p_i \to \pm\infty$. Therefore $\left\langle \frac{1}{2} a_i p_i^2 \right\rangle = \frac{1}{2}kT$, as claimed.

The condition for the proof is that the Boltzmann factor decays at the boundary, i.e. that the quadratic coefficient is positive, $a_i > 0$, and that the integration over $p_i$ can be extended to $\pm\infty$. The same manipulation applies to any variable appearing quadratically in $H$; it does not matter whether it is a coordinate or a momentum, and the potential energy $U(\{q_i\})$ plays no role in the average.

Exercise: if $H(\xi \to \pm\infty) \to \infty$, then $\left\langle \xi \frac{\partial H}{\partial \xi} \right\rangle = kT$.

Compare this with the derivation in the previous chapter. There the same identity required a careful treatment of the energy shell, integration by parts on the surface $H = E$, and the neglect of a term containing $(E-H)$. Here the shell is gone, the average is an ordinary Gaussian integral, and the proof takes three lines. This is the practical gain of the canonical ensemble.

For a Hamiltonian that is a sum of quadratic terms, $H = \sum_i A_i P_i^2 + \sum_i B_i Q_i^2$, the theorem gives $\langle H \rangle = \frac{f}{2}kT$ with $f$ the number of quadratic degrees of freedom, as derived in the previous chapter.

Ideal gases

For an ideal gas, $H = \frac{1}{2m}\sum_i p_i^2$, and

with

using

Then

and

This is the equation of state. For the internal energy, $U = A + TS = \frac{3}{2}NkT$. Eliminating $T$ in favour of $U$, i.e. $2\pi mkT = 4\pi mU/(3N)$,

Everything is the same as that from the microcanonical ensemble in the previous chapter: the equation of state, $U = \frac{3}{2}NkT$, and the entropy as a function of $V$ and $U$, including the $V/N$ under the logarithm that solves the Gibbs paradox.

Two remarks are in order. First, the factor $1/N!$ is essential; without it the entropy would contain $\log V$ instead of $\log(V/N)$, which is not extensive. Second, no factor of $h^3$ appears in the measure $\dd^3 p \dd^3 q$: a choice of the unit volume of phase space only shifts $S$ by a term proportional to $N$, i.e. by a constant per particle, and does not affect $P$, $U$, $C_V$ or any derivative of $A$. That constant cannot be fixed inside classical statistical mechanics.

The calculation is therefore classical throughout: the particles are treated as distinguishable up to permutations, and the sum over states is an integral over phase space. This describes a dilute gas at high temperature, where the particles are on average far apart compared with their quantum spread. At low temperature or high density the counting of states has to be redone with quantum statistics, and the $1/N!$ of the Gibbs factor is then replaced by the correct symmetry of the many-particle wave function.

Ising model

The Ising model is the simplest model for magnetic materials, and it exhibits a phase transition.

Fig. 3.3 shows a square lattice in two dimensions, describing a crystal. On each lattice site there is a magnet, described by $S_i = \pm 1$.

1
\begin{tikzpicture}[>=stealth, thick, font=\small, x=1.35cm, y=1.35cm]
2
  % lattice
3
  \foreach \i in {0,1,2,3} {
4
    \draw[gray!55] (0,\i) -- (3,\i);
5
    \draw[gray!55] (\i,0) -- (\i,3);
6
  }
7

8
  % spins
9
  \foreach \x/\y/\s in {
10
    0/0/1, 1/0/-1, 2/0/1, 3/0/1,
11
    0/1/-1, 1/1/-1, 2/1/1, 3/1/-1,
12
    0/2/1, 1/2/1, 2/2/-1, 3/2/-1,
13
    0/3/-1, 1/3/1, 2/3/1, 3/3/-1} {
14
      \fill (\x,\y) circle (1.6pt);
15
      \draw[->, blue] (\x,\y-0.26*\s) -- (\x,\y+0.26*\s);
16
  }
17

18
  % nearest neighbour pair
19
  \draw[red, dashed, thick] (0.15,1) -- (0.85,1);
20
  \node[red, above] at (0.5,1.12) {$\langle ij\rangle$};
21

22
  % labels
23
  \node at (1.5,3.7) {square lattice};
24
  \draw[<->] (2.55,3.05) -- (3.25,3.05);
25
  \node[above] at (2.9,3.05) {$a$};
26
  \node at (1.5,-0.75) {$S_i=\pm 1$ on every site};
27
\end{tikzpicture}

The Hamiltonian is written as

where

It is a convention to put the whole Hamiltonian in units of $kT$, because only the combination $H/kT$ enters the Boltzmann factor. With this convention the weight of a configuration is $e^{K\sum S_iS_j + h\sum S_i}$, so $K > 0$ favours parallel neighbouring spins, i.e. it is the ferromagnetic case, while $h$ favours one sign of $S_i$. The energy of a bond is $J S_i S_j = -kT K\, S_i S_j$, and the energy of a spin in the external field is $\bar h S_i = -kT h\, S_i$.

The canonical ensemble is described by $\rho \propto e^{-H/kT}$.

Omitting the Gibbs factor, the partition function is written as

Note that no $1/N!$ appears here: the lattice sites are labelled, so the configurations are distinguishable, and the sum runs over the $2^N$ configurations of the $N$ spins.

The average magnetization is

Solvable cases, with $d$ now denoting the spatial dimension:

• $d = 1$: there is no real phase transition;
• $d = 2$: too difficult, and not treated here.

The result for $h = 0$ is the following. There is a critical temperature $T_C$ such that

and near $T_C$ the magnetization behaves as

So $M$ is continuous at the transition while its derivative is not: the transition is continuous, with a jump in $\partial M/\partial T$. The exponent $\beta$ here is the usual critical exponent of the order parameter; unfortunately it is the same letter as $\beta = 1/kT$, and one has to tell them apart from the context.

How to solve, or at least approximate, such a model? The mean-field approximation replaces the interaction of a spin with its neighbours by the interaction with the average magnetization.

First, what is required of an approximation?

• After the approximation, one can perform the calculation.
• There is a clear physical meaning, and at extreme conditions, the approximation should be correct or roughly correct.

The mean-field replacement is

where $d$ is the coordination number, i.e. the number of nearest neighbours of a site; for a hypercubic lattice it equals twice the spatial dimension, so $d = 2$ for a chain and $d = 4$ for the square lattice. In the last equality the sum over pairs is counted so that each site contributes $d$ bonds, i.e. each bond is counted twice, once from each end; this is the convention that leads to $K_C = 1/d$ below.

After the replacement the Hamiltonian no longer couples different sites,

and therefore

The self-consistent condition is $M = \langle S \rangle$. For $h = 0$,

Since $\tanh x = x - x^3/3 + \cdots$, we have

The graphical solutions are shown in Fig. 3.4: the curve $\tanh(KdM)$ against the straight line $M$ itself.

1
\begin{tikzpicture}[>=stealth, thick, font=\small, x=1.15cm, y=1.15cm]
2

3
  % ---------- (a) Kd < 1 ----------
4
  \begin{scope}[shift={(0,0)}]
5
    \draw[->] (-1.65,0) -- (1.65,0) node[below] {$M$};
6
    \draw[->] (0,-1.45) -- (0,1.45);
7
    \node[above] at (0,1.45) {$\tanh(KdM)$};
8
    \draw[gray] (-1.3,-1.3) -- (1.3,1.3);
9
    \draw[blue] plot[domain=-1.5:1.5, samples=150] (\x, {1-2/(exp(0.5*\x)+1)});
10
    \fill (0,0) circle (2pt);
11
    \node at (0,-1.85) {(a) $Kd<1$: only $M=0$};
12
  \end{scope}
13

14
  % ---------- (b) Kd > 1 ----------
15
  \begin{scope}[shift={(4.8,0)}]
16
    \draw[->] (-1.65,0) -- (1.65,0) node[below] {$M$};
17
    \draw[->] (0,-1.45) -- (0,1.45);
18
    \node[above] at (0,1.45) {$\tanh(KdM)$};
19
    \draw[gray] (-1.3,-1.3) -- (1.3,1.3);
20
    \draw[blue] plot[domain=-1.5:1.5, samples=150] (\x, {1-2/(exp(4*\x)+1)});
21
    \fill (0,0) circle (2pt);
22
    \fill (0.96,0.96) circle (2pt);
23
    \fill (-0.96,-0.96) circle (2pt);
24
    \node[above left] at (0.96,0.96) {$M_0$};
25
    \node[below, yshift=-3pt] at (-0.96,-1.02) {$-M_0$};
26
    \node at (0,-1.85) {(b) $Kd>1$: $M=0$ and $M=\pm M_0$};
27
  \end{scope}
28
\end{tikzpicture}

For $Kd < 1$, there is only one solution, $M = 0$. For $Kd > 1$, there are three roots, $M = 0$ and $M = \pm M_0$, i.e. two branches.

In an equilibrium state, the free energy should reach its minimum, and one can prove that $M = \pm M_0$ is a physical solution, while $M = 0$ is only stationary. Therefore

and the phase transition takes place.

What are the features of the system around $K_C$? Let

Then

where $K_C d = 1$ has been used in going to the last line, and the term $\Delta K M^3$ has been dropped because it is of higher order. Hence

Since $K \propto 1/T$ with $J$ and $k$ fixed, $\Delta K/K_C \approx -(T-T_C)/T_C = -\tau$ to leading order near $T_C$, so $M \propto (-\tau)^{1/2}$: the mean-field value of the exponent $\beta$ is $1/2$.

Exercise: derive $\partial M/\partial \Delta K$, $U$ and $C_V$ at and around $K_C$, and discuss their behavior at and around $K_C$.

Under what conditions is the mean-field approximation almost correct?

• Higher spatial dimensions: then each spin has many neighbours, and the replacement of the neighbours by their average is less violent.
• Long-range interactions: then the mean field is not an approximation but almost the interaction itself.

It should be kept in mind that the approximation has to be justified, not assumed. In $d = 1$ the exact model has no transition at any finite temperature, and yet the mean-field condition above, with $d = 2$ for a chain, predicts $K_C = 1/2$.

Advanced methods in statistical physics:

Approximate methods

• Mean-field methods, and improved mean-field methods
• Variational methods
• Series expansions (high temperature and low temperature expansions, cumulant expansions, …)
• Renormalization group methods (real space, momentum space)

Exact methods

• Numerical solutions
• Monte Carlo simulations
• Molecular dynamics

Grand canonical ensemble

Now we allow the number of particles of the system to be non-conserved. The idea here is similar to that for deriving the canonical ensemble.

$\Gamma$ space: spanned by $(p,q,N)$, $N = 0,1,2,\dots$

The grand canonical ensemble: a collection of systems, identical in composition and macroscopic conditions, but existing in different numbers of particles and different microscopic states, without interactions with each other.

A member of the ensemble, i.e. a microscopic state of the system, is represented by a point in $\Gamma$ space. The density distribution of these points is denoted by $\rho(p,q,N)$. The arrangement is sketched in Fig. 3.1(b).

To find $\rho(p,q,N)$, we consider a canonical ensemble of $N$ particles, volume $V$, and temperature $T$. We divide the system into two subsystems,

1: $N_1, V_1$
2: $N_2 = N - N_1, V_2 = V - V_1$

and consider the special case $V_2 \gg V_1$ and $N_2 \gg N_1$.

We designate the coordinates of the $N_1$ particles in $V_1$ by $\{p_1,q_1\}$, and those in $V_2$ by $\{p_2,q_2\}$,

Note that $H(p_1,q_1,N_1)$ and $H(p_2,q_2,N_2)$ here are the same functions evaluated at different values of the arguments. This is different from the discussion of the canonical ensemble, since we now allow interchange of particles between the two subsystems: if a particle moves from $V_2$ to $V_1$, we change its coordinates from $\{p_2,q_2\}$ to $\{p_1,q_1\}$.