User:SunderB/Draft:Modelling Solids/Debye Model

Improvement from Einstein's Model: atoms are not independent - if one atom moves it causes its neighbours to move, and so on.

Collective motion in the form of sound waves - Debye wanted to quantise sound waves like Planck had quantised light

We know that the average energy of a simple harmonic oscillator is:

$$\lang E \rang = \hbar \omega[n_B(\beta \hbar \omega)+\frac{1}{2}]$$

Sound can have many differenct frequencies or 'modes'. To get total energy, we want to sum over all of these modes:

$$E_{TOT} = \sum_{modes}(\lang E_{mode}\rang) = \sum_{modes}(\hbar \omega_{mode} [n_B(\beta \hbar \omega_{mode}+\frac{1}{2})])$$

Assumptions
$$f =\frac{v}{k} \Rarr \omega=vk$$
 * Sound has three polarisation states - atoms can move in all three spatial dimensions despite the direction of travel
 * Speed of sound is independent of polarisation - inaccurate as transverse waves are usually slower than longitudinal waves, but not much more information is gained by taking this into account
 * Speed of sound is isotropic or independent of direction
 * The solid observes a linear dispersion relation:

Number of Modes in a 1D Line
$$\psi(x) = \psi_0 \sin(kx)=\psi_0 \sin(\frac{n\pi}{L}x)$$ where n is a positive integer
 * Each mode is a standing wave
 * For standing waves with nodes at each end, the wavefunction is:

Periodic Boundary Conditions
What if the line is wrapped round into a circle?

Now the boundary condition is $\psi(x) = \psi(x+L)$. Since the wave is periodic, psi can be written in the form: $$\psi(x)=e^{ikx}$$.

$$\begin{align} \psi(x) &= \psi(x+L)\\ \therefore e^{ikx} &= e^{ik(x+L)}\\ e^{ikx} &= e^{ikx}e^{ikL}\\ \therefore 1 &= e^{ikL}\\ e^{0},e^{2\pi i},e^{4\pi i}, \dots &=e^{ikL}\\ \Rarr 0, 2\pi, 4\pi, \dots &= kL\\ \Rarr k &= \frac{2n\pi}{L} \end{align} $$We want to sum up each mode

$$\therefore \lang E \rang_{TOT} = \sum_n(\lang E_n \rang)$$

since there's a unique value k for each n, we can change this to:

$$\lang E \rang_{TOT} = \sum_k(\lang E_k \rang)$$

We can change this to an integral for large values of L by considering the no. of modes in a range of k-values.

Each mode takes up $$2\pi/L$$ in 1D k-space, therefore no. of modes between $$\pm k$$ is given by:

$$\frac{k--k}{\frac{2\pi}{L}}=\frac{2kL}{2\pi}=\frac{kL}{\pi}=\frac{L}{2\pi}\int^\infty_{-\infty}dk$$

therefore, the number of modes between $$\pm\infty = \frac{L}{2\pi}\int^\infty_{-\infty}dk$$

3D Periodic Boundary Conditions
Imagine a box where if you go a distance L in any direction, you end up at a place that looks identical to where you started. We'll ignore edge effects as we're more interested in the local conditions.

In 3D, the waves can be described in the form of exponentials with vector exponents:

$$\psi(x) = e^{i\boldsymbol{k}\cdot\boldsymbol{r}}$$

where $$\boldsymbol{k} = \begin{pmatrix} k_x\\k_y\\k_z\end{pmatrix} = \frac{2\pi}{L}\begin{pmatrix} n_x\\n_y\\n_z \end{pmatrix}$$

therefore, to sum over all modes we can use:

$$\begin{align} \text{no. of modes} &= \frac{L}{2\pi}\int\Bigl(\frac{L}{2\pi}\int(\frac{L}{2\pi}\int dk_z) dk_y\Bigr) dk_x\\ \sum_k &= \Bigl(\frac{L}{2\pi}\Bigr)^3\int\int\int dk_z dk_y dk_x\\ \sum_k &= \frac{V}{(2\pi)^3}\int d \boldsymbol{k}\\ \end{align}$$

Since sound waves have 3 polarisations, we need to triple this expression:

$$\sum_{modes} = 3\sum_k = \frac{3V}{(2\pi)^3}\int d\boldsymbol{k}$$

Converting to Spherical Polar Co-ordinates
Since we're assuming speed of sound is isotropic, then we can rewrite this expression in spherical polar co-ordinates:

$$\frac{3V}{(2\pi)^3}\int d\boldsymbol{k} = \frac{3V}{(2\pi)^3}\int^{\infty}_0 4\pi k^2 dk$$

We can then use the linear velocity dispersion relation to integrate in terms of angular frequency instead of wave number:

$$\begin{align} k &= \frac{\omega}{v}\\ \therefore \frac{dk}{d\omega} &= \frac{1}{v}\\ \therefore dk &= \frac{d\omega}{v} \end{align}$$

Substituting in the expressions for k and dk, we get:

$$\begin{align} \sum_{modes} &= \frac{3V}{(2\pi)^3}\int^{\infty}_0 4\pi \frac{\omega^2}{v^3} d\omega\\ &= \frac{12V\pi}{(2\pi)^3v^3}\int^{\infty}_0 \omega^2 d\omega \end{align}$$

The integrand here is called the density of states.

In this case $$g(\omega) = \frac{12\pi V}{(2\pi)^3v^3}\omega^2$$.

If N is defined as the number of states across all frequencies in one polarisation:

$$\therefore 3N = \int^\infty_0 g(\omega) d\omega$$

We can also write the density of states in terms of a value called the Debye frequency (we'll see where this comes from later):

$$\begin{align} \omega_D &= \sqrt[3]{6\pi^2\Bigl(\frac{N}{V}\Bigr)v^3} = \sqrt[3]{6\pi^2\rho v^3}\\ \Rarr g(\omega) &= N \frac{9\omega^2}{\omega_D^3} \end{align}$$

Getting Specific Heat Capacity
Going back to our initial expression for the total average energy, we can now write it in terms of an integral of $$\omega$$:

$$\begin{align} \lang E \rang_{TOT} &= \sum_{modes}(\hbar\omega_{mode}(n_B(\beta\hbar\omega)+\frac{1}{2}))\\ \lang E \rang_{TOT} &= \int^\infty_0 \Bigl(g(\omega) \Bigl(n_B(\beta\hbar\omega)+\frac{1}{2}\Bigr)\hbar\omega\Bigr) d\omega \end{align}$$

Note that this expression is actually incorrect! The upper limit of infinity suggests that there are an infinite number of modes, and evaluating the integral will give an infinite answer, due to the term of 1/2 - this is known as the zero-point energy problem. However, for the purpose of calculating and expression for heat capacity, it doesn't have much effect as it doesn't affect the derivative. So for now, don't worry about this - we'll come back later and fix this.

Ignoring the zero-point term, we can plug in our expression for the density of states:

$$\begin{align} \lang E \rang_{TOT} &= \int^\infty_0 \Bigl(\frac{9N\omega^2}{\omega_D^3} \Bigl(n_B(\beta\hbar\omega)+\frac{1}{2}\Bigr)\hbar\omega\Bigr) d\omega\\ &= \frac{9N\hbar}{\omega_D^3} \int^\infty_0 \omega^3 n_B(\beta\hbar\omega) d\omega\\ &= \frac{9N\hbar}{\omega_D^3} \int^\infty_0 \omega^3 \Bigl( \frac{1}{e^{\beta\hbar\omega}-1} \Bigr) d\omega \end{align}$$

Substituting $$x=\beta \hbar \omega$$:

$$\begin{align} \omega&=\frac{x}{\beta\hbar}\\ \frac{d\omega}{dx}&=\frac{1}{\beta\hbar}\\ d\omega&=\frac{dx}{\beta\hbar}\\ \therefore\left\langle E\right\rangle_{\text{TOT}}&=\frac{9N\hbar}{\omega_D^3}\int_0^{\infty}\left(\frac{x}{\beta\hbar}\right)^3\left(\frac{1}{e^x-1}\right)\frac{1}{\beta\hbar}dx\\ &=\frac{9N\hbar }{\omega _D^3\beta ^4\hbar ^4}\int _0^{\infty }\frac{x^3}{e^x-1}dx\\ \end{align}$$