Skip to content
Snippets Groups Projects
Commit ca333cf8 authored by Anton Akhmerov's avatar Anton Akhmerov
Browse files

Merge branch 'move_k_space' into 'master'

Introduce k-space in lecture 2

See merge request !14
parents 339ce0a8 f5f31735
No related branches found
No related tags found
1 merge request!14Introduce k-space in lecture 2
Pipeline #14994 passed
......@@ -47,6 +47,7 @@ markdown_extensions:
- pymdownx.details
- abbr
- footnotes
- meta
extra_javascript:
- 'https://cdn.plot.ly/plotly-latest.min.js'
......
---
jupyter:
jupytext:
text_representation:
extension: .md
format_name: markdown
format_version: '1.0'
jupytext_version: 0.8.6
kernelspec:
display_name: Python 3
language: python
name: python3
---
```python
from matplotlib import pyplot
......
---
jupyter:
jupytext:
text_representation:
extension: .md
format_name: markdown
format_version: '1.0'
jupytext_version: 0.8.6
kernelspec:
display_name: Python 3
language: python
name: python3
---
```python
from matplotlib import pyplot
......@@ -20,58 +34,189 @@ _(based on chapter 2.2 of the book)_
- Write down the total energy of phonons given the temperature and the dispersion relation
- Estimate heat capacity due to phonons in high temperature and low temperature regimes of the Debye model
### Debye model
The Einstein model explained the experimental data quite well, but still slightly underestimated the observed values of $C$. Apparently the "each atom is an oscillator"-idea is too simplistic.
## Deficiency of Einstein model
The Einstein model explained the experimental data quite well, but still slightly underestimated the observed values of $C$ at low $T$. Apparently the "each atom is an oscillator"-idea is too simplistic.
We can also see that something goes wrong by comparing Einstein model to the heat capacity of silver[^1]:
```python
pyplot.rcParams['axes.titlepad'] = 20
T = np.array([1.35,2.,3.,4.,5.,6.,7.,8.,10.,12.,14.,16.,20.,28.56,36.16,47.09,55.88,65.19,74.56,83.91,103.14,124.2,144.38,166.78,190.17,205.3])
c = np.array([0.,0.,0.,0.,0.,0.,0.0719648,0.1075288,0.2100368,0.364008,0.573208,0.866088,1.648496,4.242576,7.07096,10.8784,13.47248,15.60632,17.27992,18.6188,20.33424,21.63128,22.46808,23.05384,23.47224,23.68144])
c *= 3/24.945 #24.954 is 3Nk_B
def c_einstein(T, T_E):
x = T_E / T
return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
temp = np.linspace(1, 215, 100)
fit = curve_fit(c_einstein, T, c, 500)
T_E = fit[0][0]
fig, ax = pyplot.subplots()
ax.scatter(T, c)
ax.plot(temp, c_einstein(temp, T_E), label=f'Einstein model, $T_E={T_E:.5}K$')
ax.set_ylim(bottom=0, top=3)
ax.set_xlim(0, 215)
ax.set_xlabel('$T(K)$')
ax.set_ylabel(r'$C/k_B$');
```
The einstein model works reasonably well, but it predicts a too small heat capacity at low $T$
??? question "How does $C$ predicted by the Einstein model behave at low $T$?"
When $T → 0$, $T_E/T → \infty$. Therefore neglecting $1$ in the denominator we get $C \propto \left(\frac{T_E}{T}\right)^2e^{-T_E/T}$, and the heat capacity should be exponentially small!
## Reciprocal space
Peter Debye (1884 – 1966) suggested to instead consider _normal modes_: sound waves that propagate through a solid.
> A (running) sound wave is a collective motion of atoms through this solid, where the displacement of each atom $\mathbf{\delta r}$ depends on its position $\mathbf{r}$ and time $t$ through a relation
> $$\mathbf{\delta r} = \mathbf{\delta r}_0 e^{i(\mathbf{kr}-\omega t)},$$
> with $\mathbf{\delta r}_0$ the wave amplitude, and $\mathbf{k}$ the _wave vector_ (the wave length $\lambda = 2\pi/|\mathbf{k}|$).
> Because the shape of the wave depends on time only through the factor $\exp(i\omega t)$, these waves are _normal modes_—stable oscillations of the medium.
Each normal mode has a _wave vector_ $\mathbf{k}$.
All wave vectors are points in _reciprocal space_ or _k-space_.
These waves don't have a fixed frequency $\omega_E$, but rather a _dispersion relation_
$$
\omega = v|\mathbf{k}|.
$$
Here $v$ is the *sound velocity*.
Now instead of $3N$ oscillators with the same frequency we have many oscillators with different frequencies $\omega(k) = v|\mathbf{k}|$.
This makes the total energy equal to the sum over the energies of all the oscillators:
$$
E=\sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ {\rm e}^{\hbar\omega(\mathbf{k})/{k_{\rm B}T}}-1}\right)
$$
We still have several open questions:
* Normal modes depend on the material's shape. What impact does this have on the heat capacitance?
* What $\mathbf{k}$ are possible and what aren't?
* If all $\mathbf{k}$ are possible, shouldn't $E$ be infinite?
Peter Debye (1884 – 1966) suggested to instead consider _normal modes_: sound waves that propagate through a solid with velocity $v=\omega/k$, where $k=2\pi/\lambda$ is the _wave number_. Instead of just multiplying $\bar{\varepsilon}$ by $3N$, we now use:
## Periodic boundary conditions
$$E=\int\limits_0^\infty\left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/{k_{\rm B}T}}-1}\right)g(\omega){\rm d}\omega$$
We can answer all the above questions after we observe the following.
$g(\omega)$ is the _density of states_: the number of normal modes found at each position along the $\omega$-axis. How do we calculate $g(\omega)$?
**Key idea:** $C$ shouldn't depend on material shape, and should just be proportional to its volume.
Therefore the simpler shape we consider, the easier is the calculation.
The easiest option people have invented so far is a box $L×L×L$ with **periodic boundary conditions**[^2].
#### Reciprocal space, periodic boundary conditions
Each normal mode can be described by a _wave vector_ ${\bf k}$. A wave vector represents a point in _reciprocal space_ or _k-space_. We can find $g(\omega)$ by counting the number of normal modes in k-space and then converting those to $\omega$.
This means that the displacement $\mathbf{\delta r}(\mathbf{r}$ as well as the velocity of the solid is periodic:
How far apart are the normal modes in k-space? This is determined by the boundaries of the solid. One way to treat the boundaries is by using _fixed boundary conditions_ (like a guitar string), resulting in modes $k=\pi/L$, $2\pi/L$, $3\pi/L$ etc., where $L$ is the length of the solid.
$$
\begin{align}
\mathbf{\delta r}(\mathbf{r} + (L, 0, 0)^T) &= \mathbf{\delta r}(\mathbf{r})\\
\frac{d\mathbf{\delta r}(\mathbf{r} + (L, 0, 0)^T)}{dt} &= \frac{d\mathbf{\delta r}(\mathbf{r})}{dt}\\
\end{align}
$$
In this course, however, we will exclusively use _periodic boundary conditions_, where one edge of the solid should connect seamlessly to the opposite edge. This results in:
(of course similar periodicity applies to $y$ and $z$ coordinates)
$$k=\ ..., \frac{-4\pi}{L}, \frac{-2\pi}{L}, 0, \frac{2\pi}{L}, \frac{4\pi}{L}, ...$$
Periodicity means that not all the points in $k$-space are allowed.
Instead only waves with each component $k_x, k_y, k_z$ of the $\mathbf{k}$-vector belonging to a set
$$k=…, \frac{-4\pi}{L}, \frac{-2\pi}{L}, 0, \frac{2\pi}{L}, \frac{4\pi}{L}, …$$
satisfy the periodic boundary conditions.
The three-dimensional wave vector ${\bf k}$ can be any threefold permutation of these values. All possible values of ${\bf k}$ then form a grid in k-space:
In 3D the allowed $k$-vectors form a regular grid:
![](figures/DOS_periodic.svg)
There is one allowed ${\bf k}$ per $\left(\frac{2\pi}{L}\right)^3$. The number of ${\bf k}$-values inside a sphere with radius $k$:
There is therefore exactly one allowed ${\bf k}$ per volume $\left(\frac{2\pi}{L}\right)^3$ in reciprocal space.
When we consider larger and larger box sizes $L→∞$, the volume per allowed mode becomes smaller and smaller, and eventually we obtain an integral:
$$
\sum_\mathbf{k} f(\mathbf{k}) ≈ \frac{L^3}{(2\pi)^3}\iiint\limits_{-∞}^{∞}dk_x dk_y dk_z f(\mathbf{k})
$$
## Density of states
Turning back to our problem of computing heat capacity, we get:
$$
\begin{align}
E &= \frac{L^3}{(2\pi)^3}\iiint\limits_{-∞}^{∞}dk_x dk_y dk_z × 3×\left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ {\rm e}^{\hbar\omega(\mathbf{k})/{k_{\rm B}T}}-1}\right),\\
\omega(\mathbf{k}) &= v\sqrt{k_x^2 + k_y^2 + k_z^2}.
\end{align}
$$
The factor $3$ accounts for three possible directions of displacement (wave polarizations).
To compute this integral, we observe that the integrand depends only on $|\mathbf{k}|$, and therefore switching to spherical coordinates is the way to go:
$$N=\frac{\frac{4}{3}\pi k^3}{\left(\frac{2\pi}{L}\right)^3}=\frac{Vk^3}{6\pi^2}$$
$$
\begin{align}
E &= \frac{L^3}{(2\pi)^3}\int\limits_0^{2π}d\varphi\int\limits_0^π \sin θ\;\int\limits_0^∞ k^2 dk × 3 × \left(\frac{1}{2}\hbar\omega(k)+\frac{\hbar\omega(k)}{ {\rm e}^{\hbar\omega(k)/{k_B T}}-1}\right)\\
&= \frac{L^3}{(2\pi)^3}\int_0^∞ 12 π k^2 dk \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ {\rm e}^{\hbar\omega(\mathbf{k})/{k_{\rm B}T}}-1}\right)\\
&= \frac{L^3}{(2\pi)^3}\int_0^∞ 12 π v^{-3} \omega^2 d\omega \left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/{k_{\rm B}T}}-1}\right).
\end{align}
$$
This means for the density of states $g(k)$ as a function of $k$:
In the last expression everything inside the brackets is about Bose-Einstein statistics, while all the prefactors together are specific to the problem we are studying.
$$g(k)=\frac{ {\rm d}N}{ {\rm d}k}=\frac{Vk^2}{2\pi^2}$$
We can emphasize this further by introducing a new concept, _density of states_, $g(\omega)$.
The density of states $g(\omega)$ in frequency space then becomes:
> Density of states $g(ω)$ is the number of available normal modes per infinitesimal interval $δω$.
$$g(\omega)=g(k)\frac{ {\rm d}k}{ {\rm d}\omega}=\frac{Vk^2}{2\pi^2}\frac{ {\rm d}k}{ {\rm d}\omega}$$
With this definition, our integral becomes
$$
E = ∫\limits_0^∞\left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/{k_{\rm B}T}}-1}\right)g(ω)dω,
$$
with
$$
g(ω) = \frac{L^3}{(2\pi)^3}×4 π × 3 × v^{-3} × \omega^2.
$$
We can trace back all the factors in the density of states to their origin:
In general, ${\rm d}k/{\rm d}\omega$ can be difficult to calculate; we will see more of this later. But going back to the Debye model for now, where we assume simple sound waves with $v=\omega/k$, this reduces to $g(\omega)=V\omega^2/2\pi^2v^3$. The total energy then becomes:
* $(L/2\pi)^3$ is the volume per allowed wave in the reciprocal space
* $4π$ is the area of a unit sphere, the result of integration over $d \varphi$ and $dθ$
* $ω^2$ is due to the area of this sphere being proportional to its squared radius
* $3$ is the number of possible polarizations in 3D
* $v^{-3}$ is due to $ω = v|k|$.
## Low $T$
In general, ${\rm d}k/{\rm d}\omega$ can be difficult to calculate; we will see more of this later. But going back to the Debye model for now, and using $g(\omega)=V\omega^2/2\pi^2v^3$. The total energy then becomes:
$$ E=E_{\rm Z}+\frac{3V}{2\pi^2 v_{\rm s}^3}\int\limits_0^\infty\left(\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}\right)\omega^2{\rm d}\omega$$
Here, the factor 3 comes from the fact that every wave has three polarizations (two transversal, one longitudinal). The term $E_{\rm Z}$ goes to infinity through integration. This is no problem, as it doesn't count towards the heat capacity.
Substitute $x\equiv\frac{\hbar\omega}{k_{\rm B}T}$:
$$
\Rightarrow E
= E_Z + \frac{3V}{2\pi^2 v_{\rm s}^3}\frac{\left(k_{\rm B}T\right)^4}{\hbar^3}\int\limits_0^\infty\frac{x^3}{ {\rm e}^x-1}{\rm d}x
= E_Z + \frac{3V}{2\pi^2 v_{\rm s}^3}\frac{\left(k_{\rm B}T\right)^4}{\hbar^3}\frac{\pi^4}{15}
$$
$$\Rightarrow E=E_{\rm Z}+\frac{3V}{2\pi^2 v_{\rm s}^3}\frac{\left(k_{\rm B}T\right)^4}{\hbar^3}\int\limits_0^\infty\frac{x^3}{ {\rm e}^x-1}{\rm d}x$$
Therefore we conclude that $C=\frac{ {\rm d}E}{ {\rm d}T}\propto T^3$.
The integral on the right is a constant, $\left(\frac{\pi^4}{15}\right)$ $\Rightarrow$ $C=\frac{ {\rm d}E}{ {\rm d}T}\propto T^3$.
Can we understand this without any calculation terms? Turns out we can!
#### Debye's interpolation for medium T
The above approximation works very well at low temperature. But at high temperature, $C$ should of course settle at $3k_{\rm B}$ (the Dulong-Petit value). The reason why the model breaks down, is that it assumes that there is an infinite number of harmonic oscillators up to infinite frequency.
1. At temperature $T$ only modes with $\hbar \omega \lesssim k_B T$ get thermally excited.
2. These phonons have wave vectors $|k| \lesssim k_B T /\hbar v$, and therefore their total number is proportional to the volume of a sphere with the same radius times the density of modes in $k$-space. This gives us $N_\textrm{modes} \sim (k_B T L/\hbar \omega v)^3$.
3. Each of these modes is a harmonic oscillator at a reasonably high temperature. Therefore, similarly to Einstein model, it contributes $\sim k_B$ to the heat capacity.
4. Multiplying the number of modes by the contribution of each mode we obtain $C\propto k_B (k_B T L/\hbar \omega v)^3$.
## Debye's interpolation for medium $T$
The above approximation works very well at low temperature. But at high temperature, $C$ should of course settle at $3k_{\rm B}$ (the Dulong-Petit value).
The model breaks down because it assumes that there is an infinite number of harmonic oscillators up to infinite frequency.
On the other hand we know that in total in a material with $N$ atoms, there are only $3×N$ possible normal modes.
Debye proposed an approximation: all phonons are acoustic (i.e. constant sound velocity) until a certain cut-off frequency, beyond which there are no phonons.
......@@ -87,10 +232,10 @@ $$
What determines the _Debye frequency_ $\omega_{\rm D}$?
$$
\int_0^{\omega_{\rm D}}g(\omega){\rm d}\omega=\frac{V\omega_{\rm D}^3}{2\pi^2v_{\rm s}^3}=3N,
\int_0^{\omega_{\rm D}}g(\omega){\rm d}\omega=\frac{V\omega_{\rm D}^3}{2\pi^2v_{\rm s}^3}=3N.
$$
where $N$ is the number of atom, so $3N$ the number of degrees of freedom.
(notice how $g(ω)$ comes in handy!)
Substitute in $E$, differentiate to $T$:
......@@ -101,16 +246,6 @@ $$
where $x=\frac{\hbar\omega}{k_{\rm B}T}$ and $\Theta_{\rm D}\equiv\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$, the _Debye temperature_.
```python
pyplot.rcParams['axes.titlepad'] = 20
T = np.array([1.35,2.,3.,4.,5.,6.,7.,8.,10.,12.,14.,16.,20.,28.56,36.16,47.09,55.88,65.19,74.56,83.91,103.14,124.2,144.38,166.78,190.17,205.3])
c = np.array([0.,0.,0.,0.,0.,0.,0.0719648,0.1075288,0.2100368,0.364008,0.573208,0.866088,1.648496,4.242576,7.07096,10.8784,13.47248,15.60632,17.27992,18.6188,20.33424,21.63128,22.46808,23.05384,23.47224,23.68144])
c *= 3/24.945 #24.954 is 3Nk_B
def c_einstein(T, T_E):
x = T_E / T
return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
def integrand(y):
return y**4 * np.exp(y) / (np.exp(y) - 1)**2
......@@ -119,11 +254,6 @@ def c_debye(T, T_D):
x = T / T_D
return 9 * x**3 * quad(integrand, 0, 1/x)[0]
temp = np.linspace(1, 215, 100)
fit = curve_fit(c_einstein, T, c, 500)
T_E = fit[0][0]
fit = curve_fit(c_debye, T, c, 500)
T_D = fit[0][0]
......@@ -139,6 +269,15 @@ ax.set_ylabel(r'$C/k_B$')
ax.legend(loc='lower right');
```
## Conclusions
1. Atoms in materials move in a collective fashion, forming sound waves with dispersion relation $(ω = v|\mathbf{k}|$.
2. The normal modes have a constant density of $(L/2π)^3$ in the reciprocal space.
3. The total energy and the total heat capacity are given by integrating a contribution of individual modes over the $k$-space.
4. The density of states $g(ω)$ counts the number of modes per unit frequency, and it is proportional to $ω^2$ for phonons in 3D.
5. At low temperatures phonon heat capacity is $∼ω^3$.
## Exercises
### Exercise 1: Debye model: concepts
......@@ -179,6 +318,7 @@ Explain your answer.
3. Calculate the heat capacity in the limit of high $T$ (hint: it goes to a constant).
4. At low $T$, show that $C_V=KT^{n}$. Find $n$. Express $K$ as a definite integral.
### Exercise 3: Different phonon modes
During the lecture we derived the low-temperature heat capacity assuming that all the phonons have the same sound velocity $v$.
In reality the longitudinal and transverse modes have different sound velocities (see [Wikipedia](https://en.wikipedia.org/wiki/Sound#Longitudinal_and_transverse_waves) for an illustration of different sound wave types).
......@@ -199,3 +339,7 @@ How does this change the Debye result for the heat capacity?
??? hint
Write down the total energy as an integral over $k$, then change the integration variables so that the spherical symmetry of the integrand is restored.
[^1]: Data is taken from C. Kittel, *Solid State Physics*, 2ed Wiley (1956).
[^2]: An alternative way to treat the boundaries is by using _fixed boundary conditions_ (like a guitar string), resulting in standing waves with $k=\pi/L$, $2\pi/L$, $3\pi/L$, …. This gives the same answer, but it usually more ugly, and takes more work.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment