Lo siento por la tardanza en responderte. Acabo de descubrir esta pregunta.
Preludium
En primer lugar, creo que podría haber un error o algunos no estándar de notación en su fórmula para $Q_{\mathrm{rot}}$. La rotación de los niveles de energía están dadas por
\begin{align}
E_{\mathrm{rot}} = \frac{\hbar^2}{2 I} J ( J + 1 ) \ ,
\end{align}
donde $I$ es el momento de inercia.
A continuación, la rotación constante de $B_{0} = \frac{\hbar^2}{2 I}$ (en unidades de energía) se presentó, por lo que
\begin{align}
E_{\mathrm{rot}} = B_{0} J ( J + 1 ) \ .
\end{align}
Así
\begin{align}
Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{B_0 J(J+1)}{kT}\right)} \ ,
\end{align}
mientras que su exponente contiene otro $h$ que no deberían estar ahí. Al menos esa es la notación estándar estoy acostumbrado. Para su comodidad también voy a introducir la rotación de la temperatura de $\Theta_{\mathrm{rot}} = \frac{B_{0}}{k}$, por lo que
\begin{align}
Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right)} \ .
\end{align}
Principal Respuesta
Ok, ahora que eso está resuelto, vamos a comenzar el verdadero trabajo. En realidad, es un bonito estándar de la aproximación que se utiliza aquí: el uso de la de Euler-Maclaurin fórmula
\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}
donde $f^{(k)}(a)$ $k$th derivado de la $f$ evaluado en $a$. El $b_{2i}$s son los valores absolutos de los números de Bernoulli, $b_{2} = \frac{1}{6}$, $b_{4} = \frac{1}{30}$, $b_{6} = \frac{1}{42}$, $...$.
Proporciona una manera de aproximar la suma de una función de más de una variable discreta por un análogo de la integración de más de una variable continua.
En el caso especial en la mano,
\begin{align}
f(J) = (2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right) \ ,
\end{align}
donde la suma de los límites de se $a=0$ $b=\infty$ e donde: $f(\infty) = f^{\prime}(\infty) = \ldots = 0$ esta fórmula se reduce a
\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}
El primer término es lo que se suele llamar el clásico o el de alta temperatura límite. Puede ser fácilmente calculada por hacer la sustitución $K = J(J + 1)$ $\Rightarrow$ $\mathrm{d}J = \frac{\mathrm{d}K}{2J + 1}$. Así
\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}
Ahora, sólo los derivados $f^{(k)}(0)$ necesitan ser ordenados. Aquellos que contienen términos que a continuación de la tercera orden en $\frac{\Theta_{\mathrm{rot}}}{T}$ puede calcularse como:
\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}
Hacerlo a mano es muy tedioso. Mejor usar algo como Mathematica. Conectar todo esto en el de Euler-MacLaurin fórmula y colocar todos los términos anteriores de segundo orden en $\frac{\Theta_{\mathrm{rot}}}{T}$ obtener exactamente el resultado de 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}
Como una nota al margen: Si usted tiene una asociación de molécula diatómica que usted necesita para introducir el llamado simetría número$\sigma$, lo que representa el número de distinguir las orientaciones de una molécula puede tener. Se origina a partir de la simetría limita: para las moléculas diatómicas homonucleares hay un doble eje de simetría perpendicular al eje internuclear que da lugar a dos indistinguible orientaciones, por lo que la fórmula anterior para $Q_{\mathrm{rot}}$ overcounts el número de estados cuánticos
por un factor de dos, debido a la inherente indistinguishability de la central nuclear de par. Por lo tanto, $Q_{\mathrm{rot}}$ tiene que ser dividido por $\sigma = 2$. Para moléculas diatómicas heteronucleares usted simplemente tiene $\sigma = 1$.
Con la introducción de $\sigma$ la ecuación de $Q_{\mathrm{rot}}$ se convierte en:
\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}