LECTURE 16

Systems of Interacting Particles

So far we have been considering systems of noninteracting or weakly interacting particles because these are simple to deal with. However, most systems in the real world are complicated because they have particles that interact with each other, e.g., liquids and solids. Our strategy to find the energy eigenvalues, calculate the partition function $Z$, and derive the appropriate thermodynamic quantities from $\ln Z$ still holds, but it can be difficult to find the energy eigenvalues. However there are cases of interacting systems which we can solve.

For example, at low temperatures thermodynamic properties are dominated by the low energy excited states. So we do not need to find all the energy eigenstates; we just need to describe these low energy states accurately. If the system is ordered at low temperatures, these low-lying excited states are ``collective modes.'' For example, in a crystal the atoms are located on lattice sites. The collective modes are normal modes of vibration, i.e., sound waves. Quantized sound waves are called phonons which are analogous to photons. Another example is a ferromagnet where all the spins are lined up parallel to each other in the ground state. The low energy excited states correspond to spin waves which are called magnons when quantized.

We can also sometimes solve interacting systems in the high temperature limit where $kT$ is large compared to the mean energy of interaction. In the infinite temperature limit the interactions are negligible. At high temperatures the interactions can be treated as a perturbation and one can do things like high temperature series expansions.

We will now consider some examples of interacting systems.

Lattice vibrations and normal modes

Consider a solid consisting of $N$ atoms, each of mass $m_i$ with position ${\bf r}_i=(x_{i1},x_{i2},x_{i3})$. The equilibrium position is ${\bf r}_i^{(o)}$. Each atom can vibrate about its equilibrium position. The displacement from the equilibrium position is

$$\xi_{i\alpha}\equiv x_{i\alpha}-x_{i\alpha}^{(o)}\quad\quad\quad {\rm where}\; \alpha=1,2,\;{\rm or}\;3$$ (1)

The kinetic energy of the solid is given by

$$K=\frac{1}{2}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}m_{i}\dot{x}_{... ...{2}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}m_{i}\dot{\xi}_{i\alpha}^2$$ (2)

where $\dot{x}_{i\alpha}=\dot{\xi}_{i\alpha}$ is the $\alpha$th component of the velocity of the $i$th atom.

Since the displacements are small, we can expand the potential energy $V=V({\bf r}_1,...,{\bf r}_N)$ in a Taylor series:

$$V=V_{o}+\sum_{i\alpha}\left[\frac{\partial V}{\partial x_{i\... ...partial x_{j\gamma}}\right]_{o} \xi_{i\alpha}\xi_{j\gamma}+...$$ (3)

where $i$ and $j$ go from 1 to $N$; and $\alpha$ and $\gamma$ go from 1 to 3. The derivatives are evaluated at the equilibrium positions ${\bf r}_i^{(o)}$ of the atoms. The first term $V_o$ is the potential energy when the atoms are in their equilibrium configuration. Since this is a minimum of $V$, the first derivative must vanish: $\left[\partial V/\partial x_{i\alpha}\right]_o=0$, i.e., there is no force on any atom in the equilibrium configuration. So the first term in the Taylor series vanishes. For the second term, let

$$A_{i\alpha,j\gamma}\equiv \left[\frac{\partial^{2}V}{\partial x_{i\alpha}\partial x_{j\gamma}}\right]_{o}$$ (4)

We can neglect higher order terms since the displacements $\xi$ are small. So we obtain

$$V=V_{o}+\frac{1}{2}\sum_{i\alpha,j\gamma}A_{i\alpha,j\gamma}\xi_{i\alpha}\xi_{j\gamma}$$ (5)

and the Hamiltonian becomes

$$H=V_{o}+\frac{1}{2}\sum_{i\alpha}m_{i}\dot{\xi}_{i\alpha}^2 ... ...{i\alpha,j\gamma}A_{i\alpha,j\gamma}\xi_{i\alpha}\xi_{j\gamma}$$ (6)

The kinetic energy is simple since it is just a sum of terms, each of which just involves one coordinate. But the potential energy is complicated since it involves cross terms coming from different atoms and different coordinates. This is the result of interactions. Since the potential energy is quadratic in the coordinates, we can eliminate this complication by finding the normal modes of the solid. This amounts to making a linear transformation to a new set of normal coordinates $q_s$:

$$\xi_{i\alpha}=\sum_{s=1}^{3N}B_{i\alpha,s}q_{s}$$ (7)

such that a proper choice of coefficients $B_{i\alpha,s}$ transforms the Hamiltonian to the simple (diagonal) form:

$$H=V_{o}+\frac{1}{2}\sum_{s=1}^{3N}\left(\dot{q}_{s}^{2}+\omega_{s}^2q_{s}^2\right)$$ (8)

Now we have a sum over independent harmonic oscillators with no cross terms. Each oscillator has frequency $\omega_{s}$, and its quantum mechanical energy is given by

$$\varepsilon_{s}=\left(n_s+\frac{1}{2}\right)\hbar\omega_{s}$$ (9)

where $n_s=$ 0, 1, 2 ... The total energy is the sum of these one dimensional harmonic oscillator energies:

$$E_{n_1,...,n_{3N}} = V_{o}+\sum_{s=1}^{3N}\left(n_s+\frac{1}{2}\right)\hbar\omega_{s}$$

$$= -N\eta+\sum_{s=1}^{3N}n_s\hbar\omega_{s}$$ (10)

where

$$-N\eta\equiv V_{o}+\frac{1}{2}\sum_s\hbar\omega_{s}$$ (11)

is a constant independent of $n_s$. $\hbar\omega_{s}/2$ is the zero point energy. $\eta$ represents the binding energy per atom in the solid at absolute zero.

It is now straightforward to calculate the partition function which follows what we did for the Einstein oscillators (see calculation of the Einstein specific heat in lecture 11).

$$Z = \sum_{n_1,n_2,n_3...} \exp\left(-\beta\left[-N\eta+n_1\hbar\omega_1+n_2\hbar\omega_2+...+n_{3N}\hbar\omega_{3N}\right]\right)$$

$$= e^{\beta N\eta}\left(\sum_{n_1=0}^{\infty}e^{-\beta\hbar\omega_1 ... ...right)... \left(\sum_{n_{3N}=0}^{\infty}e^{-\beta\hbar\omega_{3N}n_{3N}}\right)$$ (12)

Notice that each harmonic oscillator partition function is a geometric series. Recall that for a geometric series

$$\sum_{n=0}^{\infty}ar^n=\frac{a}{1-r}$$ (13)

So we have

$$Z=e^{\beta N\eta}\left(\frac{1}{1-e^{-\beta\hbar\omega_{1}}}\right)... \left(\frac{1}{1-e^{-\beta\hbar\omega_{3N}}}\right)$$ (14)

Now we take the logarithm to get $\ln Z$:

$$\ln Z=\beta N\eta - \sum_{s=0}^{\infty} \ln\left(1-e^{-\beta\hbar\omega_{s}}\right)$$ (15)

We can convert this sum into an integral by defining $\sigma(\omega)d\omega$ to be the number of normal modes with angular frequencies in the range between $\omega$ and $\omega+d\omega$.

$$\ln Z=\beta N\eta -\int_{0}^{\infty} \ln\left(1-e^{-\beta\hbar\omega}\right)\sigma(\omega)d\omega$$ (16)

The mean energy of the solid is

$$\overline{E}=-\frac{\partial\ln Z}{\partial\beta}=-N\eta+ \i... ...frac{\hbar\omega}{e^{\beta\hbar\omega}-1}\sigma(\omega)d\omega$$ (17)

The heat capacity at constant volume is

$$C_{V} = \left(\frac{\partial\overline{E}}{\partial T}\right)_{V} =\frac{\... ...{V} =-k_{B}\beta^{2}\left(\frac{\partial\overline{E}}{\partial\beta}\right)_{V}$$

$$= k_{B}\int_{0}^{\infty}\frac{e^{\beta\hbar\omega}} {\left(e^{\beta\hbar\omega}-1\right)^{2}}(\beta\hbar\omega)^{2}\sigma(\omega)d\omega$$ (18)

The function $\sigma(\omega)$ is determined by the normal vibrational modes of the solid. However, regardless of the exact shape of $\sigma(\omega)$, we can make some general statements about the high temperature limit. Let $\omega_{max}$ be the highest frequency of the normal mode spectrum such that:

$$\sigma(\omega)=0\quad\quad {\rm if}\;\omega > \omega_{max}$$ (19)

If the temperature is high enough such that $\beta\hbar\omega_{max}\ll 1$, then $\beta\hbar\omega\ll 1$ for all relevant $\omega$, and we can expand the exponential:

$$e^{\beta\hbar\omega}=1+\beta\hbar\omega$$ (20)

Then for $kT\gg\hbar\omega_{max}$,

$$C_{V}=k_B\int_{0}^{\infty}\sigma(\omega)d\omega=3Nk_{B}$$ (21)

since the integral over $\sigma(\omega)$ is simply the total number of modes:

$$\int_{0}^{\infty}\sigma(\omega)d\omega=3N$$ (22)

Eq. (21) is the Dulong-Petit law that we obtained earlier by applying the equipartition theorem.

Debye Approximation

The Debye approximation treats the solid as an elastic continuum and ignores the discreteness of the atoms. This is a good approximation as long as the wavelength $\lambda$ of the elastic vibration is much longer than the mean atomic lattice spacing $a$, i.e., $\lambda\gg a$. Long wavelength corresponds to low frequency modes. Let $\sigma_{c}(\omega)$ be the function describing the number of modes in an elastic continuum. Then we expect $\sigma(\omega)\approx\sigma_{c}(\omega)$ at low frequencies. This approximation will break down for $\lambda\stackrel{<}{\sim} a$.

For an elastic continuum, the dispersion relation is $\omega=c_{s}k$ where $c_{s}$ is the speed of sound. This has the same form as the dispersion relation for photons. From our previous derivation for the density of states in lecture 14, we can write

$$\sigma_{c}(\omega)d\omega=3\frac{V}{(2\pi)^3}d^3k= 3\frac{V}... ...4\pi k^2 dk\right)= 3\frac{V}{2\pi^2 c_{s}^{3}}\omega^2d\omega$$ (23)

where the factor of 3 accounts for the 3 phonon polarizations: 1 longitudinal mode and 2 modes transverse to the direction of propagation with wavevector $k$.

The Debye approximation approximates $\sigma(\omega)$ with $\sigma_D(\omega)$ defined by

$$\sigma_D(\omega)=\left\{ \begin{array}{ll} \sigma_{c}(\omeg... ...$}\\ 0 & \mbox{for $\omega > \omega_D$} \end{array} \right.$$ (24)

where the ``Debye frequency'' $\omega_D$ is chosen so that $\sigma_D(\omega)$ yields the correct total number of $3N$ normal modes:

$$\int_{0}^{\infty}\sigma_D(\omega)d\omega=\int_{0}^{\omega_D}\sigma_c(\omega)d\omega=3N$$ (25)

This normalization determines $\omega_D$:

$$\int_{0}^{\omega_D}\sigma_c(\omega)d\omega=\frac{3V}{2\pi^2c... ...D} \omega^{2}d\omega=\frac{V}{2\pi^2 c_{s}^{3}}\omega_D^{3}=3N$$ (26)

Solving for $\omega_D$ yields

$$\omega_D=c_{s}\left(6\pi^2\frac{N}{V}\right)^{1/3}$$ (27)

For typical numbers: $c_{s}\approx 5\times 10^5$ cm/sec, $a\approx(V/N)^{1/3}\approx 1\AA$, $\omega_D\approx 10^{14}$ sec $^{-1} \approx 100$ THz. The corresponding wavelength $2\pi c_{s}/\omega_{D}\sim a$. The corresponding Debye temperature $\Theta_D$ is defined by $k\Theta_D\equiv \hbar\omega_{D}$. Typically the Debye temperature is of order 300 K.

Notice that the Debye density of states is quadratic in $\omega$.

=3.0 true in

Using the Debye approximation, the heat capacity becomes

$$C_{V} = k_{B}\int_{0}^{\infty}\frac{e^{\beta\hbar\omega}} {\left(e^{\beta\hbar\omega}-1\right)^{2}}(\beta\hbar\omega)^{2}\sigma(\omega)d\omega$$

$$= k_{B}\int_{0}^{\omega_D}\frac{e^{\beta\hbar\omega}(\beta\hbar\ome... ...^{\beta\hbar\omega}-1\right)^{2}}\frac{3V}{2\pi^2c_{s}^{3}} \;\omega^{2}d\omega$$

$$= k_{B}\frac{3V}{2\pi^2\left(c_s\beta\hbar\right)^{3}} \int_{0}^{\beta\hbar\omega_D}\frac{e^x}{\left(e^x-1\right)^2}\; x^4dx$$ (28)

where $x=\beta\hbar\omega$. Using $\hbar\omega_D=k_B\Theta_D$ and

$$V=6\pi^2N\left(\frac{c_s}{\omega_D}\right)^3$$ (29)

we can write

$$C_V=9Nk_B\left(\frac{T}{\Theta_D}\right)^3\int_{0}^{\Theta_D/T} \frac{e^x}{\left(e^x-1\right)^2}\; x^4dx$$ (30)

As we have seen, at high temperatures, this reduces to the Dulong-Petit law. At low temperatures, $T\ll\Theta_D$ and we can replace the upper limit of the integral with $\infty$ and the integral becomes a constant. As a result, we can see immediately that

$$C_V\sim T^3$$ (31)

More precisely, the integral can be evaluated exactly:

$$\int_{0}^{\infty}\frac{e^x}{\left(e^x-1\right)^2}\; x^4dx=\frac{4\pi^4}{15}$$ (32)

This yields

$$C_V =\frac{12\pi^4}{5}Nk_{B}\left(\frac{T}{\Theta_D}\right)^{3}$$ (33)

Notice that $C_V\sim T^3$ at low temperatures ($T\ll\Theta_D$). This provides a reasonably good fit at low temperatures, though it may be necessary to go to temperatures as low as $T<0.02\Theta_D$. The Debye specific heat certainly gives better agreement with experiment at low temperature than the exponential temperature dependence predicted by the Einstein specific heat (see lecture 11). The Einstein specific heat makes the approximation that all the oscillators have a single frequency $\omega_E$:

$$\sigma(\omega)\approx\sigma_E(\omega)\approx 3N \delta\left(\omega-\omega_E\right)$$ (34)

Plugging this into Eq. (18) and using $\hbar\omega_E=k_B\Theta_E$ yields our previous result for the Einstein heat capacity from lecture 11:

$$C_V=3Nk_B\left(\frac{\theta_E}{T}\right)^2 \frac{e^{\theta_E/T}} {\left(e^{\theta_E/T}-1\right)^2}$$ (35)

Figure 10.2.2 in Reif shows a comparison of the Debye and Einstein specific heats. Both have an S-shape:

=3.0 true in

The Debye approximation is good for the acoustic phonon modes while the Einstein approximation is good for the high energy optical phonon modes.

=2.5 true in