# Einstein solid – numerical solution

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 3.24.

The properties of an Einstein solid can be investigated numerically for relatively small systems. The multiplicity of a solid with ${N}$ oscillators and ${q}$ energy quanta is

$\displaystyle \Omega=\binom{q+N-1}{q}=\frac{\left(q+N-1\right)!}{q!\left(N-1\right)!} \ \ \ \ \ (1)$

From this we can get the entropy

$\displaystyle S=k\ln\Omega \ \ \ \ \ (2)$

which for small systems can be worked out without using Stirling’s approximation.

The total energy of the system is

$\displaystyle U=q\epsilon \ \ \ \ \ (3)$

where ${\epsilon}$ is the energy of a single quantum, so we can get the temperature as defined from the entropy

$\displaystyle \frac{1}{T}=\frac{\partial S}{\partial U} \ \ \ \ \ (4)$

Finally, the heat capacity per oscillator is obtained from

$\displaystyle \frac{C}{N}=\frac{1}{N}\frac{\partial U}{\partial T} \ \ \ \ \ (5)$

If we produce a table of these quantities for each value of ${q}$ from 0 upwards, we can approximate the derivatives to estimate values for ${T}$ and ${C}$ using the same techniques that we used for the two-state paramagnet.

 $\displaystyle T$ $\displaystyle =$ $\displaystyle \frac{\Delta U}{\Delta S}\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\epsilon}{k}\frac{\Delta q}{\Delta\left(\ln\Omega\right)}\ \ \ \ \ (7)$ $\displaystyle \frac{C}{N}$ $\displaystyle =$ $\displaystyle \frac{1}{N}\frac{\Delta U}{\Delta T}\ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\epsilon\Delta q}{N\frac{\epsilon}{k}\Delta\left(\frac{\Delta q}{\Delta\left(\ln\Omega\right)}\right)}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{k\Delta q}{N\Delta\left(\frac{\Delta q}{\Delta\left(\ln\Omega\right)}\right)} \ \ \ \ \ (10)$

For example, we can calculate ${T\left(q\right)}$ as

$\displaystyle T\left(q\right)=\frac{U\left(q+1\right)-U\left(q-1\right)}{S\left(q+1\right)-S\left(q-1\right)} \ \ \ \ \ (11)$

and then ${\frac{C\left(q\right)}{N}}$ as

$\displaystyle \frac{C\left(q\right)}{N}=\frac{1}{N}\frac{U\left(q+1\right)-U\left(q-1\right)}{T\left(q+1\right)-T\left(q-1\right)} \ \ \ \ \ (12)$

In practice, we are actually calculating ${S/k}$, ${kT/\epsilon}$ and ${C/Nk}$.

For the case ${N=50}$ and ${q=0,\ldots,100}$, the entropy as a function of energy (actually ${U/\epsilon=q}$) is as shown

Schroeder also asks us to calculate the values for ${N=5000}$ with ${q}$ still in the same range. This dilutes the energy among a larger number of oscillators, which effectively lowers the temperature of the system since the amount of energy per oscillator is reduced. He then asks us to compare these results to the heat capacity plots in his Figure 1.14, which shows the heat capacity for one mole of lead, aluminum and diamond as functions of temperature. It’s not quite clear to me how we can compare the results of this problem with that plot, since we’re working out the heat capacity for very small numbers of oscillators (certainly much less than a mole), so the total heat capacity of our small samples is much less than the values shown in Figure 1.14. However, if the goal is to get curves that have similar shapes to those in Schroeder, we can get a crude comparison.

To do this, we need to plot the actual heat capacity and temperature rather than the dimensionless quantities we’ve used in generating the numbers. That is, we need to calculate

 $\displaystyle C$ $\displaystyle =$ $\displaystyle \frac{k\Delta q}{\Delta\left(\frac{\Delta q}{\Delta\left(\ln\Omega\right)}\right)}\ \ \ \ \ (13)$ $\displaystyle T$ $\displaystyle =$ $\displaystyle \frac{\epsilon}{k}\frac{\Delta q}{\Delta\left(\ln\Omega\right)} \ \ \ \ \ (14)$

We can vary the shape of the curve by changing the size ${\epsilon}$ of the energy quantum. For lead and aluminum, we can use the ${N=50}$ calculation to generate the following plot

The red curve uses ${\epsilon=\left(100\mbox{ K}\right)k\left(\mbox{eV K}^{-1}\right)=8.62\times10^{-3}\mbox{ eV}}$ (and extends ${q}$ out to 200), and the green curve uses ${\epsilon=\left(300\mbox{ K}\right)k\left(\mbox{eV K}^{-1}\right)=0.0259\mbox{ eV}}$ (with ${q}$ out to 100). The two curves are roughly the same shape as Schroeder’s curves for lead (red) and aluminum (green). The vertical scale for the heat capacity is, of course, much smaller than the values in Figure 1.14 since we’re dealing with only 50 particles rather than a mole.

For diamond, the heat capacity doesn’t increase much until we get to higher temperatures, so we can use the ${N=5000}$ results for that. For ${\epsilon=\left(2000\mbox{ K}\right)k\left(\mbox{eV K}^{-1}\right)=0.1724\mbox{ eV}}$ we get the following curve

which is roughly the same shape as the diamond curve in Schroeder. (The kink in the graph is due to the fact that we’re using discrete values of ${q}$.) We had to extend ${q}$ out to 200 to get this curve. Note that the vertical scale is a couple of orders of magnitude larger than the lead and aluminum plot, since we’re dealing with a hundred times as many oscillators.

I don’t know how much confidence we can place in these values for ${\epsilon}$ (probably not much), except to say that increasing ${\epsilon}$ at any given temperature tends to decrease the heat capacity. This may seem odd, since equation 13 appears to be independent of ${\epsilon}$ so it might appear that ${C}$ should not depend on ${\epsilon}$. However, if we trace the quantities in 13 back to their definitions, we see that they depend only on ${N}$ and ${q}$, so it is true that if we hold ${N}$ and ${q}$ fixed, ${C}$ does not change even if we change ${\epsilon}$. However, if we’re viewing ${C}$ as a function of temperature and not ${q}$, we note from 14 that for fixed values of ${N}$ and ${q}$, the temperature does increase as we increase ${\epsilon}$ (which makes sense, since if we have a fixed number ${q}$ of quanta and increase the energy contained in each quantum, the energy of the system increases and so does its temperature). In terms of the red/green plot above, suppose we pick a point on the red curve with heat capacity ${C_{0}}$. This point corresponds to some temperature ${T_{0}}$ which in turn is determined by values ${q_{0}}$ and ${\epsilon_{0}}$ (with everything else held fixed). If we now hold ${q_{0}}$ fixed but increase ${\epsilon}$ to a new value ${\epsilon_{1}}$, the heat capacity ${C_{0}}$ remains unchanged but the temperature increases to a new value ${T_{1}}$ which is further to the right. Thus the point on the graph moves horizontally to the right and (if we pick ${\epsilon_{1}}$ correctly) could move onto the green curve. This point at temperature ${T_{1}}$ is below the point at temperature ${T_{1}}$ on the red curve. In other words, if we want to hold the temperature constant and increase ${\epsilon}$, we must decrease ${q}$ which in turn decreases ${C}$.