1_einstein_model.md

jupyter:
  jupytext:
    text_representation:
      extension: .md
      format_name: markdown
      format_version: '1.0'
      jupytext_version: 1.0.2
  kernelspec:
    display_name: Python 3
    language: python
    name: python3

import pandas

from matplotlib import pyplot
import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import quad
import plotly.offline as py
import plotly.graph_objs as go


from common import draw_classic_axes, configure_plotting

configure_plotting()

py.init_notebook_mode(connected=True)

(based on chapter 2.1 of the book)

!!! success "Expected prior knowledge"

Before the start of this lecture, you should be able to:

- Write down the energy spectrum and partition function of a quantum harmonic oscillator
- Describe the equipartition theory
- Write down the Bose-Einstein distribution

!!! summary "Learning goals"

After this lecture you will be able to:

- Explain quantum mechanical effects on the heat capacity of solids (Einstein model)
- Compute the expected particle number, energy, and heat capacity of a quantum harmonic oscillator (a bosonic mode)
- Write down the total thermal energy of a material

Classical limit of heat capacity

Let us look at the heat capacities of different chemical elements¹:

elements = pandas.read_json('elements.json')
elements.full_name = elements.full_name.str.capitalize()
hovertext = elements.T.apply(
    lambda s: f'<sup>{s.number}</sup>{s.abbr} '
              f'[{s.full_name}{", " * bool(s.remark)}{s.remark}]'
)

fig = go.Figure(
    data=[
        go.Scatter(
            x=elements.number,
            y=elements.c / 8.314,
            mode='markers+text',
            textposition='top center',
            hovertext=hovertext,
            hoverinfo='text'
        ),
    ],
    layout=go.Layout(
        title='Heat capacity of chemical elements',
        autosize=True,
        yaxis=go.layout.YAxis(
            title='$C/k_B$',
            tick0=1,
            dtick=2,
        ),
        xaxis=go.layout.XAxis(
            title='Atomic number'
        ),
        hovermode='closest',
    ),
)

py.iplot(fig, show_link=False)

An empirical observation, also known as the law of Dulong–Petit (1819):

In most materials heat capacity per atom C \approx 3k_B

This corresponds to what we know from statistical physics and the equipartition theorem: every atom has 3 momenta and 3 coordinates. The equipartition theorem states that each of these 6 degrees of freedom contributes k_B/2 to the heat capacity.

Complication

Things start looking more complex when we study (following Einstein) the temperature dependence of the heat capacity of diamond²:

# Data from Einstein's paper
T = [222.4, 262.4, 283.7, 306.4, 331.3, 358.5, 413.0, 479.2, 520.0, 879.7, 1079.7, 1258.0]
c = [0.384, 0.578, 0.683, 0.798, 0.928, 1.069, 1.343, 1.656, 1.833, 2.671, 2.720, 2.781]

fig, ax = pyplot.subplots()
ax.scatter(T, c)
ax.set_xlabel('$T[K]$')
ax.set_ylabel('$C/k_B$')
ax.set_ylim((0, 3));

So we see that:

• At high temperature the law of Dulong–Petit works
• Strong temperature dependence of C
• At low T, C \rightarrow 0

Quantum oscillator

This can be explained by considering each atom as a quantum harmonic oscillator:

import math
from numpy.polynomial.hermite import Hermite

def ho_evec(x, n, no_states):

    """
    Calculate the wavefunction of states confined in the harmonic oscillator

    Input:
    ------
    x: numpy array of x coordinates (in units of hbar.omega)
    n: n^th bound state in the oscillator
    no_states: no of states confined

    Returns:
    --------
    Wavefunctions

    """

    # calculate hermite polynomial
    vec = [0] * no_states
    vec[n] = 1/2
    Hn = Hermite(vec)

    return ((1/np.sqrt(math.factorial(n)*2**n))*
            pow(np.pi,-1/4)*
            np.exp(-pow(x, 2)/2)*
            Hn(x))

def h0_ener(n):
    """
    Calculate the energy of nth bound state
    """
    return (n + 1/2)

x = np.linspace(-4, 4, 100) #units of hbar.omega
no_states = 7 #no of bound states confined in the quantum well

omega = 1.0 #frequency of harmonic oscillator
V = 0.5*(omega**2)*(x**2)

fig, ax = pyplot.subplots(figsize=(10, 7))

for i in range(no_states):

    ax.hlines(h0_ener(i), x[0], x[len(x)-1], linestyles='dotted', colors='k')

    ax.plot(x, ho_evec(x, i, no_states) + h0_ener(i)) #plot wavefunctions


    # annotate plot
    ax.text(x[len(x)-1], h0_ener(i)+1/4, r'$\Psi_%2i (x)$' %(i),
             horizontalalignment='center', fontsize=14)

    ax.text(1/4, h0_ener(i)+1/4, '$E_%2i$' %(i),
             horizontalalignment='center', fontsize=14)

    if i==0:
        ax.text(x[0]+1/4, h0_ener(i)/4, r'$\frac{\hbar\omega_0}{2}$',
                 horizontalalignment='center', fontsize=14)

        ax.annotate("", xy=(x[0]+1/2, h0_ener(i)-1/2),
                    xytext=(x[0]+1/2, h0_ener(i)),
                    arrowprops=dict(arrowstyle="<->"))
    elif i==1:
        ax.text(x[0]+1/4, h0_ener(i-1)+1/3, r'$\hbar\omega_0$',
                 horizontalalignment='center', fontsize=14)

        ax.annotate("", xy=(x[0]+1/2, h0_ener(i)),
                    xytext=(x[0]+1/2, h0_ener(i-1)),
                    arrowprops=dict(arrowstyle="<->"))

    ax.fill_between(x, h0_ener(i), ho_evec(x, i, no_states) + h0_ener(i), alpha=0.5)

ax.plot(x, V, 'k', linewidth=1) #plot harmonic potential

# Move left y-axis and bottim x-axis to centre, passing through (0,0)
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position(('data', 0.0))

# Eliminate upper and right axes
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

# Eliminate x and y axes labels
ax.set_yticklabels([])
ax.set_xticklabels([])

# Set x and y labels
ax.set_xlabel('X '+ r'($\sqrt{\hbar/m\omega_0}$)', fontsize=14)
ax.set_ylabel('$E/\hbar\omega_0$', fontsize=14)
ax.yaxis.set_label_coords(0.5,1)

This oscillator has an energy spectrum given by \varepsilon_n=\left(n+\frac{1}{2}\right)\hbar\omega_0 where \omega_0 is the eigenfrequency of the oscillator.

What is the thermal occupation and corresponding thermal energy stored in this oscillator? A harmonic oscillator is a bosonic mode \Rightarrow its occupation is described by the Bose-Einstein distribution: n(\omega,T)=\frac{1}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1} The Bose-Einstein distribution describes the occupation probability of a state at a given energy \hbar \omega. Using the Bose-Einstein distribution, we can calculate the expectation value of the energy stored in the oscillator
\langle E \rangle=\frac{1}{2}\hbar\omega_0+\frac{\hbar\omega_0}{ {\rm e}^{\hbar\omega_0/k_{\rm B}T}-1}

The left plot below shows the Bose-Einstein distribution vs energy. We see that low-energy states are more likely to be occupied than high-energy states. The right plot shows the increasing thermal energy in the oscillator for increasing temperature and highlights the zero-point energy \hbar\omega_0/2 that remains in the oscillator at T=0 (a consequence of the uncertainty principle). Moreover, we see that the energy in the oscillator becomes approximately constant when k_{\rm B}T\ll\hbar \omega_0. This implies that the heat capacity becomes small when k_{\rm B}T<\hbar \omega_0 and goes to zero when T\rightarrow0.

xline = [1, 1];
yline = [0, 1.5];
fig, (ax, ax2) = pyplot.subplots(ncols=2, figsize=(10, 5))
omega = np.linspace(0.1, 2)
ax.plot(omega, 1/(np.exp(omega) - 1), '-', xline, yline, 'r--')
ax.set_ylim(0, top=3)
ax.set_xlim(left=0)
ax.set_xlabel('$\hbar \omega$')
ax.set_xticks([0])
ax.set_xticklabels(['$0$'])
ax.set_ylabel('$n$')
ax.set_yticks([0,1, 2])
ax.set_yticklabels(['$0$','$1$', '$2$'])
draw_classic_axes(ax, xlabeloffset=.2)
ax.text(1.05, 0.95, r'$\hbar \omega = k_{\rm B}T$', ha='left', color='r');
temps = np.linspace(0.01, 2)
ax2.plot(temps, 1/2 + 1/(np.exp(1/temps)-1), '-', [1,1], [0, 1.1], 'r--')
ax2.set_ylim(bottom=0)
ax2.set_xlabel('$k_B T$')
ax2.set_xticks([0])
ax2.set_xticklabels(['$0$'])
ax2.set_ylabel(r"$\langle E \rangle$")
ax2.set_yticks([1/2])
ax2.set_yticklabels([r'$\hbar\omega_0/2$'])
draw_classic_axes(ax2, xlabeloffset=.15)
ax2.text(1.05, 0.65, r'$k_{\rm B}T=\hbar \omega_0$', ha='left', color='r');

We now calculate the heat capacity per atom C explicitly. To do so, we need to differentiate \langle E \rangle with respect to T. \begin{multline} C = \frac{\partial\langle E \rangle}{\partial T} = -\frac{\hbar\omega_0}{\left({\rm e}^{\hbar\omega_0/k_{\rm B}T}-1\right)^2}\frac{\partial}{\partial T}\left({\rm e}^{\hbar\omega_0/k_{\rm B}T}-1\right)\\ = \frac{\hbar^2\omega_0^2}{k_{\rm B}T^2}\frac{ {\rm e}^{\hbar\omega_0/k_{\rm B}T}}{\left({\rm e}^{\hbar\omega_0/k_{\rm B}T}-1\right)^2} =k_{\rm B}\left(\frac{\hbar\omega_0}{k_{\rm B}T}\right)^2\frac{ {\rm e}^{\hbar\omega_0/k_{\rm B}T}}{\left({\rm e}^{\hbar\omega_0/k_{\rm B}T}-1\right)^2} \end{multline} There is still one step to do, rewriting this equation more elegantly: C = k_b \left(\frac{T_E}{T}\right)^2\frac{e^{T_E/T}}{(e^{T_E/T} - 1)^2}, where we introduce the Einstein temperature T_E \equiv \hbar \omega_0 / k_B. This is a characteristic temperature when the thermal excitations "freeze out" in the harmonic oscillator. Consequently, it is also the temperature scale when the heat capacity of an Einstein solid starts significantly decreasing.

def c_einstein(T, T_E=1):
    x = T_E / T
    return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2

xline = [1, 1];
yline = [0, 1.1];
temps = np.linspace(0.01, 1.5, 500)
fig, ax = pyplot.subplots()

ax.plot(temps, c_einstein(temps)/3, '-', xline, yline, 'r--')
ax.fill_between(temps, c_einstein(temps)/3, 1, alpha=0.5)

ax.set_ylim(bottom=0, top=1.2)
ax.set_xlabel('$T$')
ax.set_ylabel(r'$C$')
ax.set_xticks([0, 1])
ax.set_xticklabels(['$0$'])
ax.set_yticks([1])
ax.set_yticklabels(['$k_B$'])
pyplot.hlines([1], 0, 1.5, linestyles='dashed')
draw_classic_axes(ax)
ax.text(1.01, 0.5, r'$T= T_E= \hbar \omega_0/k_{\rm B}$', ha='left', color='r');

The dashed line is the classical value, k_{\rm B}. Shaded area =\frac{1}{2}\hbar\omega_0, the zero point energy that cannot be removed through cooling.

This is for just one atom. In order to obtain the heat capacity of a full material, we would have to multiply C (or \langle E \rangle) by 3N, i.e. the number of harmonic oscillators according to Einstein model. The plot below shows a fit of the Einstein model to the experimental data for the heat capacity of diamond.

fit = curve_fit(c_einstein, T, c, 1000)
T_E = fit[0][0]
delta_T_E = np.sqrt(fit[1][0, 0])

fig, ax = pyplot.subplots()
ax.scatter(T, c)

temps = np.linspace(10, T[-1], 100)
ax.plot(temps, c_einstein(temps, T_E));

ax.set_xlabel('$T[K]$')
ax.set_ylabel('$C/k_B$')
ax.set_ylim((0, 3));

Conclusions

1. The law of Dulong–Petit is an observation that all materials have C≈3k_B per atom.
2. The Einstein model describes each atom in a solid as a quantum harmonic oscillator.
3. Therefore, oscillations of atoms are quantum, and they freeze out when k_BT \ll \hbar \omega_0, leaving only zero-point motion.
4. Using the Bose–Einstein distribution, we derived the thermal energy and heat capacity per atom in the Einstein model.
5. The Einstein model correctly predicts that the heat capacity drops to 0 as T\rightarrow 0.

Exercises

Quick warm-up exercises

1. Sketch the Bose Einstein distribution as a function of \omega for two different values of T
2. Sketch the heat capacity of an Einstein solid for two different values of T_E
3. What is the high-temperature heat capacity of an atom in a solid with two momentum and two spatial coordinate degrees of freedom?
4. Why is the heat capacity per atom of an ideal gas typically 3k_B/2 and not 3 k_B?
5. Explain which behaviour of the function 1/(e^{-\hbar\omega/k_BT}-1) tells you it is not the Bose Einstein distribution.

Exercise 1: Heat capacity of a classical oscillator.

Let's refresh the connection of this topic to statistical physics. You will need to look up the definition of partition function and how to use it to compute expectation values.

Consider a 1D simple harmonic oscillator with mass m and spring constant k. The Hamiltonian is given in the usual way by: H = \frac{p^2}{2m}+\frac{k}{2}x^2.

1. Compute the classical partition function using the following expression: Z = \int_{-\infty}^{\infty}dp \int_{-\infty}^{\infty} dx e^{-\beta H(p,x)}.
2. Using the solution of 1., compute the expectation value of the energy.
3. Compute the heat capacity. Check that you get the law of Dulong-Petit but with a different prefactor.
4. Explain the difference in the prefactor by considering the number of degrees of freedom.

Exercise 2: Quantum harmonic oscillator

Consider a 1D quantum harmonic oscillator. Its eigenstates are:

E_n = \hbar\omega(n+\frac{1}{2}),

1. Sketch the wave function of this harmonic oscillator for n=3.
2. Compute the quantum partition function using the following expression: Z = \sum_j e^{-\beta E_j}.
3. Using the partition function, compute the expectation value of the energy.
4. Compute the heat capacity. Check that in the high temperature limit you get the same result as in Exercise 1.3.

• What temperature can be considered high?
• What is the expectation value of n?

Exercise 3: Total heat capacity of a diatomic material

We consider a crystal of lithium, which consists of the two stable isotopes ^6Li (7.5%) and ^7Li (92.5%) in their natural abundance. Let us extend the Einstein model to take into account the different masses of these different isotopes.

1. Assume that the strength of the returning force k experienced by each atom is the same. What is the difference in the oscillation frequencies of the two different isotopes in the lithium crystal?
2. Write down the total energy stored in the vibrations of the atoms of the lithium crystal, assuming that all ^6Li atoms are in n=2 vibrational state and all ^7Li atoms are in n=4 vibrational state.
3. Write down the total energy stored in the vibrations of the atoms in the lithium crystal at a temperature T by modifying the Einstein model.
4. Compute the heat capacity of the lithium crystal as a function of T.

1. Data source: Wikipedia, mainly the CRC Handbook of Chemistry and Physics. ↩

2. The data in this plot is the same as what Einstein used, but the curve in this plot is improved compared to what Einstein did, see this blog post for the backstory. ↩