Partnership for Integration of Computation into Undergraduate Physics

# Temperature and Heat Capacity for Einstein Solids in Contact --- Previously, you built a computational model to analyze two Einstein Solids in contact and specifically to look at how the ideas thermal equilibrium and entropy can emerge from statistical combinatorics. In this activity, you will expand your previous analysis to model the temperature of two Einstein Solids in contact and the heat capacity of a single block of a few different materials. You will also compare your model of the heat capacity to empirical data. Broadly, your model should * Iterate over energy **macrostates**; Again, imagine starting with all of the quanta in solid A and then one-by-one shifting individual quanta into solid B. * Calculate the **multiplicity** and **entropy** for each macrostate (from the previous activity set) * Use macroscopic material properties (Young's Modulus) to estimate atomic properties (metalic bond strength, energy quanta) * Calculate thermodynamic quantities (**temperature** and **heat capacity**) derived from the multiplicity and the energy * Extract heat capacity data from an external file (attached to this activity) and plot those data against your model predictions. ## Exercise 1: Setting up the program If you haven't already built the model in "[Einstein Solids: Statistical Mechanics for Systems in Contact](https://www.compadre.org/PICUP/exercises/Exercise.cfm?I=354)," go do that exercise set. Call up that program to calculate the multiplicity of each macrostate for two Einstein solids---one with 150 oscillators and the other with 351 oscillators---sharing 100 quanta of energy between them. (This is already set up for you in the Jupyter code template.) **You should end up with two arrays of multiplicity values, one for each block.** Make sure that these arrays correspond to the same sequence of macrostates---that is, ```W_A[0]``` and ```W_B[0]``` should both correspond to the $q_A = 0$, $q_B=q_\text{total}$ macrostate. What percentage of the total number of oscillators are in block A? In statistical equilibrium, how many energy quanta would you expect to be in block A? (i.e. what macrostate corresponds to the greatest multiplicity? Feel free to graph this). ## Exercise 2: Calculating Temperature The statistical definition of temperature is $$ T = \dfrac{\Delta U}{\Delta S} $$ where this is regarded to be a ratio of finite differences. (This can be approximated as a derivative only in the limit of large systems with lots of energy quanta.) **We would like to calculate the temperature of each block as a function of the number of energy quanta in block A.** Given the multiplicity arrays you generated in Exercise 1, **write code that will do the following:** 1. Calculate an array of entropy values for each block. Don't forget that statistical entropy is defined as $S = k_\text{B}\ln W$ 2. Iterate through these values to create a new array of $\Delta S$ values for each block. Take a few minutes to plan how you would do this. (*Hint: Imagine iterating through each macrostate where each iteration corresponds to one block gaining an energy quanta and the other loosing a quanta. Think about your range of iteration; how many values of $S$ do you have? How many $\Delta S$'s will you have?*) 3. Create an array of $T$ values for each block. For now, set $\left|\Delta U\right| = 1$ (a single quanta), we'll change it to units of joules later in this exercise set. Don't forget to keep track of the sign of $\Delta U$ for each block! **Create a plot of $T$ vs $q_A$ for both blocks (the temperature of block B should also be plotted against the number of quanta in block A) and answer the following questions.** Don't forget that your temperature and number-of-quanta arrays will need to be of equal length when you plot them (or equal length slices of larger arrays) so pay attention to the size of your arrays. Both plots should be on the same set of axes. Include title and axes. * Where do the plots intersect? What is the physical significance of this point? Does this happen at the macrostate you would expect? * Create a second plot of $W_\text{combined}$ as a function of $q_A$. How do these plots line up? ## Exercise 3: Calculating $\Delta U$ The Einstein Solid models the chemical bonds between atoms as "spring-like" quantum harmonic oscillators. https://www.glowscript.org/#/user/matterandinteractions/folder/matterandinteractions/program/12-wells-oscillator The energy of a single quantum is: $$ \begin{align} \Delta U &= 1\text{ quantum}\\ &= \hbar \omega\\ &= \hbar \sqrt{\dfrac{k_\text{s}}{m}}\\ &= \hbar \sqrt{\dfrac{4k_\text{s,bond}}{m}} \end{align} $$ $m$, in this case is the mass of a single atom and $k_\text{s,bond}$ is the effective "spring stiffness" of the chemical bond between neighboring atoms. The factor of 4 comes from the fact that in each dimension, each atom has a bond on both sides (for a factor of 2) and each spring-like bond is "cut-in-half" with the other half servicing the neighboring atom (a second factor of two)---that is, atoms in the bulk of a solid will have a larger oscillation frequency than a simple diatomic molecule would have with an equivalent bond. (We've already taken into account the fact that each atom has 3 oscillators by setting $N = 3N_\text{atoms}$ in the combinatoric equation) But how can we calculate $k_\text{s,bond}$? It turns out that the **macroscopically observable modulus of linear elasticity ([Young's Modulus](https://en.wikipedia.org/wiki/Young's_modulus))** is a macroscopic manifestation of microscopic atomic bond strength. $$ Y = \dfrac{\text{stress}}{\text{strain}} = \dfrac{(F/A)}{(\Delta L/L)} $$ Rearranging to recreate the form of Hooke's law: $$ \begin{align} F &= \left(\dfrac{YA}{L}\right)\Delta L\\ &= k_\text{s}\Delta L \end{align} $$ This is a macroscopic expression, but if we imagine zooming into the stretching of a single bond, then we can write $$ \begin{align} F &= \left(\dfrac{YA_\text{atom}}{L_\text{bond}}\right)\Delta L_\text{bond}\\ &= k_\text{s,bond}\Delta L_\text{bond} \end{align} $$ So, one way to estimate the effective bond stiffness is to just multiply young's modulus by the cross-sectional area of the volume occupied by a single atom and then divide by the atomic diameter (or just simply multiply by the atomic diameter). **Write code to do the following:** * Use the values given below to calculate the bond stiffness and size of an energy quanta for both lead and aluminum * Edit your calculation of temperature to feature these new calculations of $\Delta U$. You should end up with an array of temperature values for each material. * For each material, plot the temperature of block A as a function of $q_A$. These should all be in the same plot. Be sure to include proper labeling, etc. | - | Aluminum | Lead | |---|----|----| |Young's Modulus | $68.9\times10^9$ N/m$^2$| $16\times10^9$ N/m$^2$| |Density| 2700 kg/m$^3$| 11340 kg/m$^3$ | |Atomic Mass | $4.48\times10^{-26}$ kg/atom| $3.44\times10^{-25}$ kg/atom| **Answer the following questions** * Do these material values comport with your knowledge that lead is heavy but "soft" while aluminum is comparatively light but hard? * Does the equilibrium temperature depend on the material of the blocks? (Assume that both are the same material; You can try making one lead and one aluminum to see what happens but since these materials have different sized quanta, then the model's exchange of individual quanta would not conserve energy.) ## Exercise 4: Empeirical Heat Capacity Data In the next few parts, we'll calculate the heat capacity of both an aluminum block and a lead block and compare these models to real data. From here on out, **focus just on block A**, which we will take to be alternatively aluminum and lead. **Download the csv file of heat capacity data attached to this exercise set and save this file in the same folder as your program.** If you cannot access this csv file, a copy is publicly available on GitHub https://github.com/QuantumBrandon/Statistical-Mechanics-Public/blob/master/SpecificHeatData.csv. **Write code to have your program open this file and extract the data.** If you do not know how to do this, this is a good opportunity to find documentation online. Plot the heat capacity vs. temperature. Make this plot look nice. ## Exercise 5: Numerical Prediction of C The definition of heat capacity is $$ C = \dfrac{\Delta U}{\Delta T} $$ where, as with temperature, this is taken to be a ratio of discrete differences. **Write code to calculate $C(T)$ for both Aluminum and Lead. You should end up with an array of C values for each material. Plot these predictions as functions of temperature alongside the empirical data.** Both sets of plots should be on the same set of axes. Include title, axes, and a legend. As in your calculation of $T$, don't forget to keep track of your iteration range and that the ``plot`` function requires you to input arrays of equal length (or equal length slices of larger arrays). Some things to keep in mind: * Look back at your solution to Exercise 2 if you need hints on how to proceed. * You only need to consider one block gaining energy quanta, but repeat this for each material. * Your numerical predictions will have dimensions of heat capacity per ($N_\text{oscillators}/3$) atoms while the empirical data has dimensions of heat capacity per mole; make sure that the dimensions of your model are consistent with the data before plotting. (If you get lost in the unit conversion, pay attention to the y-axis). * Your "bond spring stiffness" is an estimation based on the material average; feel free to "tune" your prediction in order to fit the empirical data. * If you would like to test this model against other materials, csv files attached to this set contain data for diamond and silicon. **Answer the following questions** * How well does the Einstein Model of a solid reproduce the empirical data? How does this compare to the classical prediction ($C = 3N_\text{A}k_\text{B}$ for all monatomic solids, independent of temperature)? What does this imply about the validity of the model? What are some aspects of the empirical data that don't seem to be captured by your model? What are some aspects of the physical situation that your model doesn't capture (at least as far as heat capacity goes)? ## Exercise 6: Analytic Prediction of C You may have noticed that the numerical plot doesn't extend down close to zero; this is partly a limitation of keeping our values for $N$ and $q$ small enough that the combinatoric functions don't overload the processor. As an alternative, we can turn to the analytic solution for heat capacity of the Einstein solid: $$ C(T) = 3N_\text{A}k_\text{B}\dfrac{e^{\hbar\omega/k_\text{B}T}}{\left(e^{\hbar\omega/k_\text{B}T} - 1\right)^2}\left(\dfrac{\hbar\omega}{k_\text{B}T}\right)^2 $$ **Write code that does the following:** 1. Write a function that takes as input $\omega$ and $T$ and outputs a value for the heat capacity. 2. Create an array of temperature values between 1 and 400 and input these into this heat capacity function along with the oscillation frequency for aluminum. 3. Plot the analytic solution against both the empirical data for aluminum and the numerical model for aluminum. This plot should include appropriate labels, etc. Comment on these. What are some aspects of the empirical data that don't seem to be captures by your model?