Return to Answer

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}

Main Answer

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{B_0 J(J+1)}{kT}\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}

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}

Main Answer

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}

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}

Main Answer

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{B_0 J(J+1)}{kT}\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}

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}

Main Answer

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{B_0 J(J+1)}{kT}\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}

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}

Main Answer

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_{i}}{(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_{i}$s are the Bernoulli numbers, $b_{1} = \frac{1}{6}$, $b_{2} = \frac{1}{30}$, $b_{3} = \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{B_0 J(J+1)}{kT}\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_{i}}{(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}

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}

Main Answer

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{B_0 J(J+1)}{kT}\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}