We aim to calculate the specific heat of a diatomic molecule as a function of temperature. We know the atomic energy levels due to rotation: e j(j+1), j=0,1,2,... so the partition function, Z, provides an exact answer in theory. Two problems are soon displayed: mathematics does not know the exact answer to the infinite sum for Z AND if the sum is calculated numerically, the terms in the sum span a range too large to be handled by your PC's version of "real" numbers (a problem you also saw in calculating Omega, which was solved by writing special purpose code). Also note that if the infinite sum is terminated to a finite number of terms, at high enough temperature, the Boltzmann factors all become one---meaning the states are equally occupied---so further increase in temperature results in no change in average energy so the specific heat approaches zero. SO no fixed sum cut-off will work at all temperatures. The solution we adopt below is to use more terms in the sum at higher temperatures, and finally use the integral approximation for the sum at "high" temperatures. (The integral approximation does not well-approximate the contributions of small j states, but at high temperatures the occupancy of those states doesn't much change with temperature and hence does not contribute much to the specific heat.) We select as our sample diatomic molecule H2 (molecular hydrogen). I note that H2 actually becomes a liquid at 20 K, so specific heats near or below 20 K are not accurate. Additionally real H2 has additional complexities due to the fermion/boson issues briefly mentioned in class. So this model, which will require a lot of effort, is not the final story. (But for us it will suffice.) At low temperatures, the large j states (with energy e j(j+1)) will have very small Boltzmann factors, and hence will not much contribute to the sum. For example j=5 at T=100 K has Boltzmann: Exp[-87.6 5*6/100] 3.86137 10^-12 where for e we have used e=87.6*k, and, as usual, k is the Boltzmann constant (which ended up cancelling with kT)...see below: In general e is calculated using the moment of inertia, I: e=hbar^2/(2 I^2) the nuclei of H2 are separated by d=0.7416e-10 m I=2 mp*(d/2)^2 * let I=2*mp*(0.7416e-10/2)^2 * ? hbar^2/(2*I) .1208968783469029E-20 * ? res/kb 87.56528297796090 We are interested in e/k=87.6 K as that is what frequently appears in Z. As usual we try an leave the value of variables unfixed until numerical results are required as in a plot. Note we initially use b=beta=1/kT, as derivatives wrt b define the average energy, but then switch to T as derivatives wrt T define specific heat. The specific heat we calculate is that divided by N k, as that is what we've seen is often an integer or half-integer. For "low" temperatures we quit the sum at j=4. z1=Sum[(2 j +1) Exp[-e b j(j+1)],{j,0,4}] avg1=-D[Log[z1], b] avg1=avg1 /. b->1/(k T) c1=D[avg1,T] c1 = c1 /. {e-> 87.6 k} c1=Simplify[c1/k] Plot[c1,{T,10,100}] Some remarks about this code: the notation "/." means substitute, for example we replace b with 1/(k T) and we put in the value for e appropriate for H2. We wanted to scale out from c the N k....we never bothered to put in the N (we've been calcualting the average energy of a molecule), the k is removed by Simplify[c1/k]. Note that this completely removes k; we never have to enter its numerical value. Howework Questions: Explain what/why of the following bit of code: avg1=-D[Log[z1], b] c1=D[avg1,T] Try the following plots: Plot[c1,{T,1,10}] Plot[c1,{T,100,1000}] The low temperature plot fails because of the smallest possible machine real number; the high temperature plot shows a declining specific heat because of the finite sum cut-off. For temperatures above 100 K, we need more terms: z2=Sum[(2 j +1) Exp[-e b j(j+1)],{j,0,8}] avg2=-D[Log[z2], b] avg2=avg2 /. b->1/(k T) c2=D[avg2,T] c2 = c2 /. {e-> 87.6 k} c2=Simplify[c2/k] Plot[c2,{T,100,250}] Note that in the range 80-200 c1 & c2 closely match: Plot[{c1,c2},{T,80,200}] c2 has "small number" issues below T=60 K Note that c2 also goes wrong a high temperatures: Plot[c2,{T,250,1000}] For the highest temperatures we'll use the integral: z3=Integrate[Exp[-e b u],{u,0,Infinity},Assumptions->{e b >0}] avg3=-D[Log[z3], b] avg3=avg3 /. b->1/(k T) c3=D[avg3,T] c3 = c3 /. {e-> 87.6 k} c3=Simplify[c3/k] Note that c3 equals one (as was shown in the text). We can define a function that selects one of these results as appropriate for the temperature allc[T_]:=If[T<100,c1, If[T<300,c2,c3]] Note the unusual way the independent variable is denoted when defining a function: underscore: T_ If I wanted to define a hypotenuse: hyp[x_,y_]:=Sqrt[x^2+y^2] Note2: there is a difference between ":=" and "=", but it is of no current concern. You would then use the function without the underscore (i.e., you only use it when defining it). Plot[allc[T],{T,10,200},PlotRange->All] Turn in hardopy of this plot! Note that this plot matches Figure T5.2: the specific heat goes to zero at low temperatures ("freeze out"), goes to one at high temperatures (equipartition theorem), and over-shoots by a bit in between.