Assignment 4: Heat capacity, chemistry, and 1D materials

Exercise 1 - Debye says hi

a. Explain what are the differences between the Einstein model and the Debye model, explicitly stating the assumptions of the Debye model and what issues they sought to correct

The Debye model considers only sound waves in a material (i.e. one assumes a linear dispersion $\omega = v k$) which oscillate up to a maximum frequency $\omega_{\mathrm{cutoff}} = \omega_d$, the Debye frequency, which exists to ensure that the total number of vibrational modes in the system is correct (equal to the number of degrees of freedom in the system).

b. Show that the "cutoff" frequency introduced by Debye, $\omega_{cutoff}$ is indeed the Debye frequency

The cutoff frequency is defined by

$$\int_0^{\omega_{cutoff}} g(\omega)d\omega = 3N$$

We also have the relation the

$$\sum_{modes}\rightarrow \frac{V}{(2\pi)^3}\int d\mathbf{k}$$

So if we set out to count the modes:

$$\begin{align*} N_{modes} = 3 \times \sum_{k\,modes} & = 3 \times \frac{V}{(2\pi)^3}\int d\mathbf{k} \\ & = 3 \times \frac{V}{(2\pi)^3} 4\pi \int_0^{k_D} k^2 dk \\ & = 3 \frac{4\pi V}{(2\pi)^3} \frac{1}{v_s^3}\int_0^{\omega_D} \omega^2 d\omega \\ & = \int_0^{\omega_D}\frac{12\pi V}{(2\pi)^3 v_s^3} \omega^2 \,d\omega \equiv \int_0^{\omega_{D}} g(\omega)d\omega \end{align*}$$

and Debye demands the $N_{modes} = 3N$, so then we have the two expressions

$$\int_0^{\omega_{cutoff}} g(\omega)d\omega = 3N = \int_0^{\omega_{D}} g(\omega)d\omega$$

and clearly this means $\omega_{cutoff}=\omega_D$

c. Write an expression for the number of modes in a two-dimensional system, and thus determine the \emph{Debye wavenumber} (the wavenumber which corresponds to the Debye frequency).

In two dimensions, there should be $2N$ modes (we implicitly assume that the different sound polarisations propagate with the same velocity, and thus

$$\begin{equation} 2N = 2 \frac{A}{2\pi^2} \int_0^{|k| = k_d} \mathrm{d}^2 k \end{equation}$$

where $A$ is the area of our periodic system ($L \times L$). The integral is simply the area of a circle with a radius $k_d$, so it evaluates to

$$\begin{equation} 2N = \frac{2A}{(2\pi)^2}(\pi k_d^2) \end{equation}$$

and ultimately obtains

$$\begin{equation} k_d = \sqrt{4\pi n} \end{equation}$$

with $n = N/A$ (the two-dimensional density).

d. Provide a brief discussion of which elements you would expect to have a high Debye temperature, and which elements you would expect to have a low Debye temperature.

Note: your reasoning is important, not the actual elements!

The Debye model considers only sound waves in a material (i.e. one assumes a linear dispersion $\omega = v k$) which oscillate up to a maximum frequency $\omega_{\mathrm{cutoff}} = \omega_d$, the Debye frequency, which exists to ensure that the total number of vibrational modes in the system is correct (equal to the number of degrees of freedom in the system).

Exercise 2 - Sommerfeld

a. In your own words, explain what is the Fermi energy, Fermi temperature and the Fermi surface

The Fermi energy is the chemical potential at $T=0$, the Fermi temperature is the associated temperature $T_F = E_F/k_{\textrm{B}}$, and is the temperature at which thermal energy would dictate the properties of a material to the same degree as the fact they are fermions. The Fermi surface is the surface in momentum space which separates the filled and unfilled states at zero temperature.

b. Write an expression for the number of states for a gas of free electrons in three dimensions and use this to calculate the Fermi wavenumber and Fermi Energy

$$N = \frac{2V}{(2\pi)^3} \int_{k<k_F} \mathrm{d}\mathbf{k}$$

The factor of two is very important. The integral is just the volume of a sphere of radius $k_F$, this is, $4\pi k_F^3/3$ which then yields

$$k_F = (3\pi^2 N/V)^{1/3}$$

and thus the Fermi energy is

$$E_F = \frac{\hbar^2(3\pi^2 N/V)^{2/3}}{2m}$$

c. Using the previous result, estimate the value of the Fermi energy for Caesium

Hint: You will need to research values for the properties of Caesium

Caesium has a density of $1.93\textrm{g/cm^3}$ and a mass of approximately $133 \textrm{amu}$ (132.9055) and thus the density of atoms is

$$N/V = (1.93 \textrm{g/cm^3}) \times (10^2 \textrm{cm/m})^3 \times (\frac{1}{133}\textrm{mol/g}) \times (N_A \textrm{atoms/mol}) = 8.7 \textrm{atoms/m^3}$$

and with one valence electron, $N/V$ for electrons is identical and we can plug it into the expression for $E_F$ from part (ii) which returns $E_F = 2.48 \times 10^{-19} \textrm {J} = 1.6 \textrm{eV}$

d. Now that you have done this for a three-dimensional system, let's consider the two-dimensional system: find an expression for the density of states $g(\varepsilon)$ at the Fermi surface of a $2D$ free-electron gas, and relate it to the Fermi energy.

Hint: Repeating the steps of the three-dimensional case above would be a good start.

For a two-dimensional gas

$$\begin{equation} N = 2A \int_{k<k_F} \frac{\mathrm{d}\mathbf{k}}{(2\pi)^2} = \frac{2A}{(2\pi)^2}\pi k_F^2 \end{equation}$$

where (again) $A$ is the area and the integral is over a circle with radius $k_F$ (c.f. exercise 1 b., only a factor of two different due to polarisation). Therefore

$$\begin{equation} k_F = (2\pi N/A)^{1/2} \end{equation}$$

and the Fermi energy is

$$\begin{equation} E_F = \frac{\hbar^2 (2\pi N/A)}{2m} = \frac{\hbar^2\pi N/A}{m} \end{equation}$$

The density of states $g(\varepsilon)$ can be acquired in a few ways, but the simplest is to use the relation

$$\begin{equation} N = \frac{A}{2\pi} k^2 \end{equation}$$

for $k<k_F$, and we know the dispersion relationship is

$$E=\frac{\hbar^2k^2}{2m}$$

and thus

$$\begin{align*} g(\varepsilon) = \frac{dN}{d\varepsilon} & = \frac{d}{d\varepsilon} \frac{A}{2\pi} \frac{2mE}{\hbar^2} \\ & = \frac{A m}{\pi \hbar^2} = N/E_F \end{align*}$$

which is surprisingly, independent of energy.

Exercise 3 - Bonding: not LCAO

a. In you own words, explain why ionic bonds occur, and what properties one would expect from and ionic solid.

Ionic bonds occur as an electron is transferred from one atom to another, and the resulting ions attract each other. Typical properties are due to the nature of the bond being very strong, with materials having high melting temperatures, and usually being hard, brittle, and electrically insulating. They are also mostly soluble, but this is not of so much relevance here.

b. The (first) ionisation energy of sodium is roughly $5.14~\textrm{eV}$, and the electron affinity of chlorine is roughly $3.62~\textrm{eV}$, and the bond length between the two atoms when a sodium chloride molecule is formed is roughly $0.236~\textrm{nm}$. Assuming that all of the cohesive energy is due to the Coulomb interaction, calculate the bonding energy.

The cohesive energy is related to the bond distance $d$ via

$$\begin{equation} E_{coh} = \frac{e^2}{4\pi\epsilon_0 d}. \end{equation}$$

Using the value $d= 2.36 \textrm{Å}$ one finds a cohesive energy of $6.10~\textrm{eV}$ and thus the total binding energy is

$$\begin{equation} E = -5.14 + 3.62 + 6.10 = 4.58~\textrm{eV} \end{equation}$$

c. The measured value of the bonding energy of sodium chloride is $4.26~\textrm{eV}$. How does this compare to your value above? Provide an explanation for any discrepancies.

The calculated value above is slightly larger than the measured value with the reason for the discrepancy being there must be a repulsive fore between the ions (otherwise the would collapse into a singularity!) in addition to the Coulomb attraction, therefore reducing the size of the cohesive (binding) energy.

d. In our discussion of bonding, we did not explicitly discuss van der Waals bonds. Research what is the nature of the van der Waals bond, explicitly describing the origin of the attractive force formation and reason as to why the force is of the form $R^{-7}$ What properties might one expect for a system bonded by van der Walls bonds?

Van der Waals forces are from correlated dipole fluctuations. If the electron is a given fixed position, there is a dipole moment $\mathbf{p} = e\mathbf{r}$ where $\mathbf{r}$ is the vector from the electron to the proton. With the electron "orbiting" (i.e. in an eigenstate), the average dipole moment is zero. However, if an electric field is applied to the atom, the atom will develop a polarisation (i.e., it will be more likely for the electron to be found on one side of the nucleus than on the other). We write

$$\begin{equation} \mathbf{p} = \chi \mathbf{E} \end{equation}$$

with $\chi$ the polarisability. Now, let us suppose we have two such atoms, separated a distance $r$ in the $\hat{x}$ direction. Suppose one atom momentarily has a dipole moment $\mathbf{p}_1$ (for definiteness, suppose this dipole moment is in the $\hat{z}$ direction). Then the second atom will feel an electric field

$$\begin{equation} E = \frac{p_1}{4\pi\epsilon_0 r^3} \end{equation}$$

in the negative $\hat{z}$ direction. The second atom then, due to its polarisability, develops a dipole moment $p_2 = \chi E$ which in turn is attracted to the first atom. The potential energy between these two dipoles is

$$\begin{equation} U=\frac{-\left|p_{1}\right|\left|p_{2}\right|}{4 \pi \epsilon_{0} r^{3}}=\frac{-p_{1} \chi E}{\left(4 \pi \epsilon_{0}\right) r^{3}}=\frac{-\left|p_{1}\right|^{2} \chi}{\left(4 \pi \epsilon_{0} r^{3}\right)^{2}} \end{equation}$$

corresponding to a force which is attractive and proportional to $1/r^7$. \hfill\mrk{1} Note that while for a single isolated atom $\langle p \rangle = 0$ the result is proportional instead to $\langle |p|^2 \rangle \sim \langle |r|^2 \rangle$ with $r$ the position of an electron, is nonzero.

Exercise 4 - One-dimesnional oscillations

a. Explain what is meant by a normal mode and a phonon

People tend to be confused by phonons, so to explicitly connect the two and explain why phonons are bosons: each classical normal mode of vibration corresponds to a quantum mode of vibration which can be excited multiple times. A single mode may be occupied by a single phonon, or it may be occupied with multiple phonons corresponding to a larger amplitude oscillation. The fact that the same state may be multiply occupied by phonons means that phonons must be bosons.

b. Show that the dispersion relation for longitudinal oscillations of an infinite one-dimensional chain of identical atoms, assuming mass $m$, spring constant $\kappa$, and lattice spacing $a$ is given by

$$\omega = 2\sqrt{\frac{\kappa}{m}}\left|\sin\left(\frac{ka}{2}\right)\right|$$

The equation of motion of the $n^{\textrm{th}}$ particle along the chain is given by

$$\begin{equation} m \ddot{x}_{n}=\kappa\left(x_{n+1}-x_{n}\right)+\kappa\left(x_{n-1}-x_{n}\right)=\kappa\left(x_{n+1}+x_{n-1}-2 x_{n}\right) \end{equation}$$

where $na$ is the equilibrium position of the $n^{\textrm{th}}$ particle. Looking for solutions of the form

$$\begin{equation} x_n = A e^{i\omega t - i k n a} \end{equation}$$

one obtains

$$\begin{equation} \begin{aligned} -\omega^{2} m e^{i \omega t-i k n a} &=\kappa e^{i \omega t}\left(e^{i k(n+1) a}+e^{i k(n-1) a}-2 e^{i k n}\right) \\ \omega^{2} m &=2 \kappa(\cos (k a)-1) \end{aligned} \end{equation}$$

which can be recast as

$$\begin{equation} \omega = \sqrt{\frac{2\kappa}{m}\left(\cos(ka)-1\right)} = 2\sqrt{\frac{\kappa}{m}}\left|\sin\left(\frac{ka}{2}\right)\right| \end{equation}$$

c. Show that a the mode with wavevector $k$ is equivalent to the mode $k + 2\pi/a$

Nothing difficult about this one:
$$\begin{equation} e^{-i(k+2\pi/a)na} = e^{-i k n a} e^{-i 2\pi n} = e^{-i k n a} \end{equation}$$

d. Assuming periodic boundary conditions, what is the spacing of adjacent modes, and how many different modes are there?

If one assumes periodic boundary conditions, then $k = 2 \pi m/L$, but $k$ is identified with $k + 2 \pi /a$ so that there are therefore exactly $N = L/a$ different normal modes.

e. Show that the phase and group velocities are $$v_p = \frac{2\sqrt{\frac{\kappa}{m}}\left|\sin\left( \frac{ka}{2} \right)\right|}{k}$$

and $$\begin{align*} v_g & = a\sqrt{\frac{\kappa}{m}} \cos\left(\frac{|k|a}{2}\right) \\ &=\frac{a \omega_0}{2}\sqrt{1-\frac{\omega^2}{\omega_0^2}} \end{align*}$$ where $\omega_0 = 2\sqrt{\kappa/m}$. It is expected that you will demonstrate the validity of both expressions for $v_g$.

The phase velocity is calculated via

$$v_p = \omega/k = \frac{2\sqrt{\frac{\kappa}{m}}\left|\sin\left( \frac{ka}{2} \right)\right|}{k}$$

and the group velocity is

$$v_g = \frac{d\omega}{dk} = a\sqrt{\frac{\kappa}{m}} \cos\left(\frac{|k|a}{2}\right).$$

We can then set $\omega_0 = 2\sqrt{\kappa/m}$ which gives

$$\frac{a}{2} \omega_0 \left\rvert \cos\left( \frac{ka}{2} \right) \right\rvert = \frac{a}{2} \omega_0 \sqrt{ 1−\sin^2\left(\frac{ka}{2}\right)}$$

and it follows directly that

$$\begin{equation} v_g = \frac{a \omega_0}{2}\sqrt{1-\frac{\omega^2}{\omega_0^2}} \end{equation}$$

f. Show that the density of states $g(\omega)$ can be written as

$$g(\omega) = \frac{2 N}{2 \pi \sqrt{(\kappa / m)-(\omega(k) / 2)^{2}}}$$

Hint: find an expression for $g(k)$ and then use the chain rule to relate $dN/dk$ to $dN/d\omega$. Results from part e. may also be of some use.

In the normal way, this model has baked into it the idea of a periodic system, so this means that we will have discretised $k$ states, which will be equally spaced and thus The density of states is uniform in $k$. Specifically, if there are $N$ sites in the system, there are $N$ modes of oscillation, and therefore the number of states per unit length is $g(k) = dN/dk = L/(2\pi) = Na/(2\pi)$ where $L$ is the length of the system. \mrk{1} We then have

$$\begin{align*} g(\omega) &=d N / d \omega=(d N / d k)(d k / d \omega)=\frac{N a}{2 \pi} \frac{1}{v_{\text {group }}} \\ &=\frac{N}{2 \pi} \frac{1}{\sqrt{\kappa / m} \cos (|k| a / 2)} \\ &=\frac{2 N}{2 \pi \sqrt{(\kappa / m)-(\omega(k) / 2)^{2}}} \end{align*}$$

where we have used

$$\left(\frac{\omega}{2}\right)^2 + \left(\frac{v_g}{a}\right)^2 = \frac{\kappa}{m}$$

which follows directly from our final expression in part e.

Exercise 5 - Using a computer to generate heat

a. In the previous question, you may have shown that the density of states for an 1D harmonic chain is given by

$$g(\omega) = \frac{2 N}{2 \pi \sqrt{(\kappa / m)-(\omega(k) / 2)^{2}}}$$

Using Python, produce a plot of $g(\omega)$ for meaningful values of $\omega$, and comment on the the physical significance of this plot. Explicitly, discuss why the density of states would be greater at one frequency as compared to another with reference to the harmonic chain and its quantum analogue.

Plotting the density of states

N = 1
fig, ax = plt.subplots()
w = np.linspace(0, np.pi/a - 1.15, 300);
g = (2*N)/(2 * np.pi * np.sqrt((kappa/m)-(w/2)**2));
ax.plot(w, g, 'C0');
ax.set_xlabel(r'$\omega$');
ax.set_ylabel('$g(\omega) [N/(\pi \sqrt{k/m})]$');
ax.set_title('Density of states')
plt.yticks([N/(np.pi * np.sqrt(kappa/m)), 2, 3], [1, 2, 3]);
ax.set_ylim(0,3)
plt.tight_layout();
plt.show()

The density of states tells you the number of states per unit frequency or energy. The reason that the density of states increases at the Brillouin zone boundary is that oscillation modes are essentially standing waves, meaning that the group velocity of the wave (see part e from question 4) vanishes. So for low $k$, a wave will freely propagate so we end up without many states per unit energy, but as the velocity slows down, we have many more states per unit energy. A common analogy is cars: on a highway with a large velocity, there are few cars per unit length, but when you geat near a boundary - which would correspond to a traffic jam - there are many more cars per unit length. Not perfect, but it is sufficiently illustrative. In terms of the quantumn analogue, I am not sufficiently clear-headed about what I wanted here: in general, the translation of the harmonic chain to the quantum analogue is that oscillators are quantised which gives rise to phonons in the system, but the above arguments hold here.

b. Using $g(\omega)$, produce a plot of the heat capacity versus temperature for the 1D harmonic chain. Does this plot match your expectations? Ensure to justify your answer in physical terms.

The energy stored in the chain is given by

$$\begin{equation} U = \int d\omega~g(\omega)~\hbar \omega \left(n_{\textrm{B}}(\omega)+1/2\right) \end{equation}$$

and so the heat capacity is $\partial U/\partial T$. The factor of 1/2 can be ignored as it is a temperature independent constant and thus will vanish upon differentiation.

Plotting this requires numerical integration of the above equation, and one should get a plot similar to that as shown below, with the values on the $x$ axis changing with the parameters $\kappa$, $m$, and $a$.

Computing the heat capacity

a = 1e-10
def integrand(w):
   b = 1/T
   nb = 1/(np.exp(b * w) - 1)
   g = (2*N)/(2 * np.pi * np.sqrt((kappa/m)-(w/2)**2))
   return g * w * nb
temp = np.linspace(0.01, 5, 150)
U = []
for T in temp:
   U.append(integrate.quad(integrand, 0, 2*np.sqrt(kappa/m))[0])
dt = temp[1]-temp[0]
dudt = np.gradient(U, dt)
fig, ax = plt.subplots()
ax.plot(temp, dudt)
ax.set_xlabel('$T$ [K]')
ax.set_ylabel('$C/k_\mathrm{B}$')
ax.set_title('Specific heat')
plt.show()

The plot should look like something that satisfies all the content we have covered in class, namely, it will be a constant ($k_\mathrm{B}$) for large $T$ and as $T$ decreases it will decay exponentially.