# How to make approximation of rotational partition function of diatomic linear molecules?

Using the rigid rotor approximation to the level energies, and such other appropriate assumptions, we can approximate rotational partition function, $$Q_{\mathrm{rot}}$$, of linear molecules as follows: $$\begin{eqnarray} Q_{\mathrm{rot}} & = & \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{E_J}{kT}\right)}\\ & \simeq & \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{hB_0J(J+1)}{kT}\right)} \\ \end{eqnarray}$$ In the book, Microwave Molecular Spectra (Gordy & Cook, 1984), Chapter 3, Equation 3.64, the authors made approximation as: $$Q_{rot}\simeq\frac{kT}{hB_0}+\frac{1}{3}+\frac{1}{15}\left(\frac{hB_0}{kT}\right)+\frac{4}{315}\left(\frac{hB_0}{kT}\right)^2$$ Also, there is another approximation made by McDowell, R. S. 1988, Journal of Chemical Physics, 88, 356: $$Q_{rot}\simeq\frac{kT}{hB_0}\exp\left(\frac{hB_0}{3kT}\right)$$ And we can see from the picture that $$Q_{\mathrm{rot}}\simeq\frac{kT}{hB_0}+\frac{1}{3}$$ is also a good approximation. Then I can guess how to derive McDowell's approximation, but I cannot derive Gordy & Cook's approximation.

What I have got are as follows:

Since $$n\exp(-x(n+1)n)+\int_{(n+1)n}^{(m+1)m}{\exp(-xy)}\,\mathrm dy+(m+1)\exp(-x(m+1)m) \\ \leq\sum_{J=n}^{m}{(2J+1)\exp(-x(J+1)J)} \\ \leq\sum_{i=n^2}^{m^2+2m}{\exp(-xi)}$$ hence, $$\int_{0}^{\infty}{\exp(-xy)}\,\mathrm dy=\frac{1}{x} \\ \leq\sum_{J=0}^{\infty}{(2J+1)\exp(-x(J+1)J)} \\ \leq\sum_{i=0}^{\infty}{\exp(-xi)}=\frac{1}{1-\exp(-x)}$$

i.e., $$\frac{kT}{hB_0}

While the latter one approximates $$\frac{1}{1-\exp\left(-\frac{hB_0}{kT}\right)}\approx \frac{kT}{hB_0}+\frac{1}{2}+\frac{1}{12}\frac{hB_0}{kT}-\frac{1}{720}\left(\frac{hB_0}{kT}\right)^3$$ when $$T$$ is rather high.

Is there anyone knowing how to derive Gordy & Cook's approximation? Please tell me and I appreciate it very much.

Besides, another little question, are such approximations still useful as nowadays computers are greatly powerful?

• This question appears to be off-topic because it is better suited for Physics.SE Jul 19 '14 at 12:49
• I only just discovered the question so I couldn't answer it earlier. I've provided an answer which hopefully explains everything. Jul 19 '14 at 18:30
• @Martin, still thank you and the answer has come:) Jul 20 '14 at 2:52

Sorry for the late answer. I just discovered this question.

## Preludium

Firstly, I think you might have an error or some non-standard notation in your formula for $Q_{\mathrm{rot}}$. The rotational energy levels are given by

\begin{align} E_{\mathrm{rot}} = \frac{\hbar^2}{2 I} J ( J + 1 ) \ , \end{align}

where $I$ is the moment of inertia. Then the rotational constant $B_{0} = \frac{\hbar^2}{2 I}$ (in units of energy) is introduced, so

\begin{align} E_{\mathrm{rot}} = B_{0} J ( J + 1 ) \ . \end{align}

Thus

\begin{align} Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{B_0 J(J+1)}{kT}\right)} \ , \end{align}

whereas your exponent contains another $h$ which should not be there. At least that's the standard notation I am used to. For convenience I'll also introduce the rotational temperature $\Theta_{\mathrm{rot}} = \frac{B_{0}}{k}$, so that

\begin{align} Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right)} \ . \end{align}

Ok, now that that's sorted out, let's start the real work. It is actually a pretty standard approximation that is used here: They use the Euler-Maclaurin formula

\begin{align} \sum_{n = a}^{b} f(n) &= \int_{a}^{b} f(n) \, \mathrm{d} n + \frac{1}{2} \bigl( f(b) + f(a) \bigr) + \sum_{i=1}^{\infty} (-1)^{i} \frac{b_{2i}}{(2i)!} \Bigl( f^{(2i - 1)}(a) - f^{(2i - 1)}(b) \Bigr) \end{align}

where $f^{(k)}(a)$ is the $k$th derivative of $f$ evaluated at $a$. The $b_{2i}$s are the absolute values of the even Bernoulli numbers, $b_{2} = \frac{1}{6}$, $b_{4} = \frac{1}{30}$, $b_{6} = \frac{1}{42}$, $...$. It provides a way to approximate the summation of a function over a discrete variable by an analogous integration over a continuous variable.

In the special case at hand,

\begin{align} f(J) = (2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right) \ , \end{align}

where the summation limits are $a=0$ and $b=\infty$ and where $f(\infty) = f^{\prime}(\infty) = \ldots = 0$ this formula boils down to

\begin{align} Q_{\mathrm{rot}} = \sum_{J = 0}^{\infty} f(J) &= \int_{0}^{\infty} f(J) \, \mathrm{d} J + \frac{1}{2} f(0) + \sum_{i=1}^{\infty} (-1)^{i} \frac{b_{2i}}{(2i)!} f^{(2i - 1)}(0) \end{align}

The first term is what is usually called the classical or the high-temperature limit. It can be easily calculated by making the substitution $K = J(J + 1)$ $\Rightarrow$ $\mathrm{d}J = \frac{\mathrm{d}K}{2J + 1}$. So

\begin{align} \int_{0}^{\infty} f(J) \, \mathrm{d} J &= \int_{0}^{\infty} (2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right) \, \mathrm{d} J \\ &= \int_{0}^{\infty} \exp\left(-\frac{\Theta_{\mathrm{rot}} K}{T}\right) \, \mathrm{d} K \\ &= \frac{T}{\Theta_{\mathrm{rot}}} \end{align}

Now, only the derivatives $f^{(k)}(0)$ need to be sorted out. Those that contain terms below the third order in $\frac{\Theta_{\mathrm{rot}}}{T}$ can be calculated to be:

\begin{align} f(0) &= 1 \\ f^{(1)}(0) &= 2 - \frac{\Theta_{\mathrm{rot}}}{T} \\ f^{(3)}(0) &= - 12 \frac{\Theta_{\mathrm{rot}}}{T} + 12 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \\ f^{(5)}(0) &= 120 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - 180 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + 30 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{5} \ . \end{align}

Doing this by hand is very tedious. Best you use something like Mathematica. Plugging all this into the Euler-MacLaurin formula and dropping all terms above second order in $\frac{\Theta_{\mathrm{rot}}}{T}$ you get exactly the result of Gordy & Cook

\begin{align} Q_{\mathrm{rot}} &\approx \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{2} \cdot 1 + (-1) \frac{1}{6 \cdot 2!} \left( 2 - \frac{\Theta_{\mathrm{rot}}}{T} \right) \\ &\quad + (-1)^{2} \frac{1}{30 \cdot 4!} \left( - 12 \frac{\Theta_{\mathrm{rot}}}{T} + 12 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \\ &\quad + (-1)^{3} \frac{1}{42 \cdot 6!} \left(120 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - 180 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + 30 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{5} \right) \\ &= \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{2} - \frac{1}{6} + \frac{1}{12} \frac{\Theta_{\mathrm{rot}}}{T} - \frac{1}{60} \frac{\Theta_{\mathrm{rot}}}{T} \\ &\quad + \frac{1}{60} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \frac{1}{252} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \\ &= \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{3} + \frac{1}{15} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \end{align}

As a side note: If you have a homonuclear diatomic molecule you need to introduce the so called symmetry number $\sigma$ which represents the number of indistinguishable orientations a molecule can have. It originates from symmetry constrains: for homonuclear diatomic molecules there is a two-fold axis of symmetry perpendicular to the internuclear axis which gives rise to two indistinguishable orientations, so the above formula for $Q_{\mathrm{rot}}$ overcounts the number of available quantum states by a factor of two because of the inherent indistinguishability of the nuclear pair. Thus, $Q_{\mathrm{rot}}$ needs to be divided by $\sigma = 2$. For heteronuclear diatomic molecules you simply have $\sigma = 1$. With the introduction of $\sigma$ the equation for $Q_{\mathrm{rot}}$ becomes:

\begin{align} Q_{\mathrm{rot}} &\approx \frac{1}{\sigma} \left( \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{3} + \frac{1}{15} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \right) \\ &= \frac{T}{\sigma \Theta_{\mathrm{rot}}} \left(1 + \frac{1}{3} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{1}{15} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} \right) \right) \end{align}

• +1, this is a really exemplary answer. I just wanted to provide a link to an article that I found useful when studying the connection between the Euler-Maclaurin formula and Bernoulli numbers in the context of studying the calculus of finite differences. Jul 19 '14 at 20:32
• @GregE. Thanks for the kind comment and the link. Looks interesting. Jul 19 '14 at 20:41
• @Philipp, thank you. I really appreciate it. As for $h$, is it due to c.g.s units? C.g.s units are still commonly used in astronomy, and it might still need a long time for the conversion to S.I. units. Jul 20 '14 at 2:46
• @puresky Glad to hear that my answer is helpful. Unfortunately, I'm not familiar enough with c.g.s units to tell you whether they are the reason for the additional $h$. Might well be but I can't say for sure. This would be a question better suited for Physics.SE I think. By the way, I added some information about homonuclear diatomic molecules and the symmetry number that needs to be introduced when you want to calculate the partition sum for them at the end of my answer. Jul 20 '14 at 12:20