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 ipywidgets import interact, FloatSlider, Layout

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 prerequisites"

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 theorem
- Write down the Bose-Einstein distribution

!!! summary "Learning goals"

After this lecture you will be able to:

- Explain the effect of the quantum harmonic oscillator on the heat capacity of solids (the Einstein model)
- Compute the expected occupation number, energy, and heat capacity of a quantum harmonic oscillator (a bosonic mode)
- Write down the total internal energy of an Einstein solid

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 various 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 the heat capacity per atom C \approx 3k_B

This corresponds to what we know from statistical physics. Under the assumption that the atomic potential is parabolic, the equipartition theorem states that each degree of freedom contributes k_B/2 to the heat capacity. Because we consider a 3D solid, each atom contains 3 spatial and 3 momentum degree of freedom. Therefore, the total heat capacity per atom is C = 3k_B.

Complication

However, at low temperatures we see that the heat capacity of diamond drops below the prediction of the law of Dulong–Petit. This suggests that we need a different model to describe the heat capacity at low temperatures. Following Einstein's reasoning² we will now explain this puzzle.

# 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))
ax.set_title('Heat capacity of diamond as a function of the temperature');

We observe that:

• We obtain the law of Dulong–Petit at high temperatures
• C depends strongly on the temperature
• As T \rightarrow 0, C \rightarrow 0

the Einstein model

The equipartition theorem assumed that each atom can be modeled as a classic harmonic oscillator. However, at low temperatures this led to a discrepancy in the heat capacity between the law of Dulong-Petit and the observed heat capacity. To explain the behaviour of C at low temperatures, Einstein treated each atom as a quantum harmonic oscilator, instead of a classic harmonic oscillator. Einstein also assumed that each atom oscillates with the same frequency \omega_0. To summarize, the Einstein model of a solid is characterised by the following two assumptions:

• The atoms are independent quantum harmonic oscillators
• Each atom has the same frequency \omega_0

For simplicity, we consider a 1D quantum harmonic oscillator (it turns out that the heat capacity in the 3D case is simply 3 times the heat capacity of the 1D case).

The energy levels and wavefunctions of a 1D quantum harmonic oscillator are shown below.

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)

The energy spectrum of a quantum harmonic oscillator is given by

E_n=\left(n+\frac{1}{2}\right)\hbar\omega_0

where \omega_0 is the eigenfrequency of the oscillator, which is the same for each atom.

What is the thermal occupation and corresponding thermal energy stored in this oscillator? A harmonic oscillator is a bosonic mode, therefore its occupation number is described by the Bose-Einstein distribution: n_{BE}(\omega,T) = \frac{1}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1} The Bose-Einstein distribution describes the expected occupation number 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
\begin{align} \langle E \rangle &= \frac{1}{2}\hbar\omega_0+\hbar\omega_0 \langle n \rangle \\ &= \frac{1}{2}\hbar\omega_0+\frac{\hbar\omega_0}{{\rm e}^{\hbar\omega_0/k_{\rm B}T}-1} \end{align} Here we have substituted the Bose–Einstein distribution n_{BE}(\omega,T) for the expected occupation number \langle n \rangle.

The left plot below shows the Bose-Einstein distribution as a function of the energy \hbar omega. We observe that low-energy states have a higher expected occupation number than high-energy states.

The right plot shows the expected value of the energy as a function of the temperature. We see that for high T, \langle E \rangle is linear in T. For low T, there is insufficient energy to excite the atoms into an excited state and therefore all atoms occupy the ground state (n = 0). Hence, \langle E \rangle converges towards a constant value of \hbar omega_0/2, which is the zero-point energy. 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\ll\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), '-')
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_{BE}(\omega,T)$')
ax.set_yticks([0,1, 2])
ax.set_yticklabels(['$0$','$1$', '$2$'])
ax.set_title(r'$n_{BE}$ as a function of $\hbar \omega$')
# 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), '-', [0.55,0.55], [0, 0.7], '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$'])
ax2.set_title(r'$\langle E \rangle$ as a function of $T$')
draw_classic_axes(ax2, xlabeloffset=.15)
ax2.text(0.65, 0.35, r'$k_{\rm B}T=\hbar \omega_0$', ha='left', color='r');

Having found an expression for \langle E \rangle as a function of T, we can now calculate the heat capacity per atom C explicitly. To do so, we need to differentiate \langle E \rangle with respect to T. \begin{align} 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{align} We can rewrite this equation into C = k_b \left(\frac{T_E}{T}\right)^2\frac{e^{T_E/T}}{(e^{T_E/T} - 1)^2}, where we have introduced the Einstein temperature T_E \equiv \hbar \omega_0 / k_B.

This is the characteristic temperature for which the thermal excitations "freeze out" in the harmonic oscillator. This means that there is not enough thermal energy to excite the harmonic oscillators into an excited state, which leaves them in the ground state. Consequently, it is also the temperature scale for which the heat capacity of an Einstein solid starts significantly decreasing. We also observe that the Einstein temperature depends on the eigenfrequency \omega_0. As different materials have a different eigenfrequency \omega_0, the Einstein temperature is a material-dependent parameter.

# Defining variables
y_line = [0, 0.92];
x_max = 2.5
temps = np.linspace(0.01, x_max, 750)

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

def einstein_temp_plot(T_E):
    # Creating plot
    fig, ax = pyplot.subplots(figsize = (12, 9))
    pyplot.hlines([1], 0, x_max, linestyles = 'dashed', label = r'Classical')
    ax.plot(temps, c_einstein(temps, T_E)/3, '-', label = r'Einstein model')
    ax.plot([T_E, T_E], y_line, 'r--')
    ax.fill_between(temps, c_einstein(temps, T_E)/3, 1, alpha = 0.5)
    ax.set_title('Heat capacity of the Einstein model and the equipartition theorem')
    ax.set_ylim(bottom = 0, top = 1.2)
    ax.set_xlabel('$T$')
    ax.set_ylabel(r'$C$')
    ax.set_xticks([0])
    ax.set_xlim((0, x_max))
    ax.set_xticklabels(['$0$'])
    ax.set_yticks([1])
    ax.set_yticklabels(['$k_B$'])
    ax.text(T_E+0.05, 0.5, r'$T= T_E= \hbar \omega_0/k_{\rm B}$', ha='left', color='r')
    ax.legend(loc = 'lower right')
    pyplot.show()
    return()

# Creating the slider
layout_slider = Layout(width = '725px', height = '20px')
float_slider = FloatSlider(
    value = 1, 
    min = 0.1, max = 2, step = 0.05,
    description = r'$T_E$',
    continuous_update = True,
    orientation = 'horizontal',
    readout_format = '.1f',
    layout = layout_slider)

# Creating the plot
control = interact(einstein_temp_plot, T_E = float_slider)

The horizontal dashed line is the classical value, k_{\rm B}. The shaded area is the difference between the classical value k_B and the value predicted by the Einstein model. Integrating over the shaded area yields \frac{1}{2}\hbar\omega_0, which is the zero-point energy of the oscillator, which cannot be extracted from the system. The vertical dashed line depicts the Einstein temperature T_E, at which the heat capacity C \approx 0.92 k_B.

The calculations above were done for a single harmonic oscillator. In order to obtain the heat capacity of a full material, we would have to multiply C (or \langle E \rangle) by 3N. The N comes from the number of harmonic oscilators in the solid and the factor 3 originates from the dimensionality of the system. ??? 1D versus 3D harmonic oscillator

    One 3D (quantum) harmonic oscillator is equivalent to three independent 1D (quantum) harmonic oscillators. Therefore, we can simply multiply by a factor of 3.

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, label = r'Experimental value')

temps = np.linspace(10, T[-1], 100)
ax.plot(temps, c_einstein(temps, T_E), label = r'Einstein model');
ax.set_title(r'Emperical and predicted heat capacity of diamond as a function of $T$')
ax.set_xlabel('$T[K]$')
ax.set_ylabel('$C/k_B$')
ax.set_ylim((0, 3));

Although the Einstein model fits the experimental data quite well, it still deviates from the experimental data in the low temperature regime.

Conclusions

1. The law of Dulong–Petit is an observation that all materials have C \approx 3k_B per atom.
2. The Einstein model describes each atom in a solid as an independent quantum harmonic oscillator with the same eigenfrequency \omega_0.
3. Using the Bose–Einstein distribution, we derived an expression for \langle E \rangle and C as a function of the temperature.
4. At sufficiently low T, the thermal excitations freeze out, resulting in \langle E \rangle = \hbar \omega_0/2.
5. The Einstein model correctly predicts that the heat capacity drops to 0 as T\rightarrow 0.

Exercises

Quick warm-up exercises

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

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)}. where \beta = 1/k_B T
2. Using the solution of 1., compute the expectation value of the energy.
3. Calculate the heat capacity. Does it depend on the temperature?

Exercise 2: Quantum harmonic oscillator

Consider a 1D quantum harmonic oscillator. Its energy eigenvalues are:

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

1. Compute the partition function using the following expression: Z = \sum_j e^{-\beta E_j}.
2. Using the partition function found in 2.1, compute the expected value of the energy.
3. Compute the heat capacity. Check that in the high temperature limit you get the same result as in Exercise 1.3.
4. Plot the found heat capacity and roughly indicate in the plot where the Einstein temperature is.
5. What is the expected value of n?

Exercise 3: Total heat capacity of a diatomic material

One of the assumptions of the Einstein model states that every atom in a solid oscillates with the same frequency \omega_0. However, if the solid contains different types of atoms, it is unreasonable to assume that the atoms oscillate with the same frequency. One example of such a solid is a lithium crystal, 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 mode and all ^7Li atoms are in n=4 vibrational mode.
3. In the case where the osccilators can occupy any vibrational mode, 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. ↩