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.