please see the attached file.

We need to evaluate the partition function for the rotational degrees of freedoms for one molecule in the high temperature limit, Z_{rot}.

We can write the total partition function as:

Z = Z_{rot}*Z_{trans}

where Z_{trans} is the partition function for the translational degrees of freedom. Because the thermodynamical quantities are functions of (derivatives of) logarithms of Z, they can be written as a sum of a purely translational part and a rotational part.
We can write any of the thermodynamic quantity X as:

X = X_{rot} + X_{trans}

The translational part is just:

U_{trans} = 3/2 NkT

S_{trans} = Sackur-Tetrode formula

F_{trans} = U_{trans} - T S_{trans}

H_{trans} = U_{trans} + PV = 5/2 NkT

G_{trans} = F_{trans} + PV = F_{trans} + NkT

The chemical potential mu follows from the fact that G = N mu (let me know if you don't know how this is derived)

So, all we have to do is to calculate the quantities X_{rot} for all the above quantities.

The eigenvalues for the operator L^2 are given by h-bar^2 l(l+1). The energy is:

L^2/(2I) and this can thus take the values:

h-bar^2/(2I) l (l+1)

The degeneracy of the energy levels is (2l+1). If we put h-bar^2/(2I) equal to epsilon, which is given as 0.00025 eV we can write the partition function (for one molecule) as:

Z = sum over l from zero to infinity (2l+1) exp[-beta epsilon l(l+1)].

For high temperatures the term beta epsilon will change only slowly and we can approximate the summation by a integral. A systematic way to expand the summation in powers of beta epsilon is as follows:

One can write a summation:

Sum from n1 to n2 of f(n) as:

Integral over x from n1 to n2 of f(x)dx + 1/2 [f(n1) + f(n2)] +

sum over r from 1 to infinity B_{2r}/(2r)! [f^(2r-1)(n2) - f^(2r-1)(n1)].

This is the Euler-McLaurin summation formula, which I explain in detail in the attached PDF file (I give a heuristic derivation), see formula 10 . The numbers B_{2r} are the Bernoulli numbers. For n2--> infinity and n1 = 0, the first few terms of the above expansions gives you:

Integral over x from 0 to infinity of f(x)dx + 1/2 f(0) -1/12 f'(0) + 1/720 f'''(0)

I our case:

f(l) = (2l+1) exp[-beta epsilon l(l+1)]

The first term is the integral from zero to infinity, which yields:

1/(beta epsilon)

The second term yields:

1/2 f(0) = 1/2

The third term yields:

-1/12 f'(0) = -1/12[2-beta epsilon]

The fourth term yields (to order beta epsilon):

1/720 f'''(0) = 1/720 (-12 beta epsilon)

And you can easily see that the later terms won't contribute to order beta epsilon.

If you sum everything up you get the result:

Z_{rot} = 1/(beta epsilon) + 1/3 + beta epsilon/15

You now use the standard formulae to obtain the thermodynamic quantities from Z_{rot} (remember that this is per molecule, so we must multiply by N):

U_{rot} = -Nd Log(Z)/d beta = N[1/beta - epsilon/3 - beta epsilon^2/45]

F_{rot} = -N/beta Log[Z_{rot}] = 1/beta Log[beta epsilon] - epsilon/3 - beta epsilon^2/90

These expressions are all accurate to order beta epsilon^2, the neglected term is of order beta^2 epsilon^3.

From F_{rot} and U_{rot} you solve for S_{rot}:

F_{rot} = U_{rot} - T S_{rot}

Z_{rot} does not depend on the volume so there is no contribution to the pressure. This means that the enthalpy is just:

H_{rot} = U_{rot} + PV = U_{rot} + NkT

And also:

G_{rot} = F_{rot} + NkT

and the contribution of the rotational degrees of freedom to the chemical potential follows from:

G_{rot} = N mu_{rot}.