Ho-kay... I see how this is done. It's going to involve layers upon layers of small terms expansions, so I'm not going to derive the whole series shown because I don't want to keep track of that many terms or do that many integrals. So, here goes:
The first thing to do is, as you guess, convert the sum into an integral:
\ln Z_\pi = -\nu \sum_{\mathbf{k}} \ln\left(1 - e^{-\beta \varepsilon_k}\right) \rightarrow -\frac{\nu \Omega}{(2\pi)^3}\int_0^{2\pi} d\phi \int_{-1}^{1}d (\cos\theta) \int_0^{\infty} dk~k^2 \ln\left(1 - e^{-\beta \varepsilon(k)}\right)
Now, do the angular integrals since the energy depends only on the magnitude, and integrate the magnitude integral by parts:
-\frac{4\pi\nu\Omega}{(2\pi)^3}\left\{ \left.\frac{k^3}{3}\ln\left(1-e^{-\beta\sqrt{k^2+m_\pi^2}}\right)\right|^{\infty}_0 + \frac{\beta}{3}\int_0^{\infty}dk k^4 \frac{e^{-\beta \sqrt{k^2+m_\pi^2}}/\sqrt{k^2 + m_\pi^2}}{1-e^{-\beta \sqrt{k^2-m_\pi^2}}}\right\}
The integrated out piece from the integral vanishes, leaving you with a nasty, nasty integral that would probably cause Mathematica to punch you in the face if you tried to enter it into it. So, we do the next best thing, because it's the only thing we know how to do: assume m_\pi is small compared to k and expand as a series. To start, I'll rewrite the integral a little bit:
\ln Z_\pi \simeq \frac{\nu \Omega \beta}{6\pi^2}\int_0^{\infty}dk~k^3\frac{\left(1+ \left(\frac{m_\pi}{k}\right)^2\right)^{-1/2}}{e^{\beta k\left(1 + \left(\frac{m_\pi}{k}\right)^2\right)^{1/2}} - 1}
So, we're going to be expanding the square root terms. I assume you're familar with this so I'm not going to do those steps to save myself some latexing:
\frac{\nu \Omega \beta}{6\pi^2}\int_0^{\infty}dk~k^3\frac{\left(1 - \frac{1}{2}\left(\frac{m_\pi}{k}\right)^2 + \cdots\right)}{e^{\beta k}e^{\beta k\left( \frac{1}{2}\left(\frac{m_\pi}{k}\right)^2 + \cdots \right)} - 1}
Now we're going to want to expand the second expontential in the denominator there, as its argument is small (at least if we add on the additional assumption that \beta m_\pi is small.
\frac{\nu \Omega \beta}{6\pi^2}\int_0^{\infty}dk~k^3\frac{\left(1 - \frac{1}{2}\left(\frac{m_\pi}{k}\right)^2 + \cdots\right)}{e^{\beta k}\left[1 + \left(\beta k\left( \frac{1}{2}\left(\frac{m_\pi}{k}\right)^2 + \cdots \right)\right) + \frac{1}{2}\left(\beta k\left( \frac{1}{2}\left(\frac{m_\pi}{k}\right)^2 + \cdots \right)\right)^2 + \cdots \right] - 1}
So, as you can see, this thing is getting uuuuugggglyyy. So I'm just going to go up to O(m_\pi^2). Splitting the numerator into two terms, we see that in the second term I can neglect the expansion in the numerator as that'll give me terms of order higher than I really want. In the first term I need to expand the denominator. Factor out an exp(\beta k) - 1 from the denominator and expand.
The zeroth order term is easy to evaluate:
\frac{\nu \Omega \beta}{6\pi^2}\int_0^{\infty}dk~\frac{k^3}{e^{\beta k} - 1} = \frac{\nu \Omega \beta}{6\pi^2}\beta^{-4} 3! \zeta(4) = (\nu \Omega \beta)T^4 \frac{\pi^2}{90}
The integral is evaluated noting that it is proportional to the polylogarithm \mbox{Li}_{4} evaulated at 1, which is equivalent to the Reimann zeta function evaluated at 4. the proportionality factor is 3! in this case. (The factor of beta^(-4) came from defining x = beta*k to de-dimensionalize the integral). Evidently Boltzmann's constant is set to 1 in this book.
Now for the second order terms:
-\frac{m_\pi}{2}\int_0^{\infty}dk~k\left[\frac{e^{\beta k}\beta k}{\left(e^{\beta k} -1\right)^2} + \frac{1}{e^{\beta k} - 1} \right]
The first term came from expanding the denominator of the ugly part above. The second was from the second order term of the numerator. Cleaning up the integral a bit gives the total second order contribution as
-\frac{\nu \Omega \beta T^4}{12 \pi^2}\frac{m_\pi^2}{T^2}\int_0^{\infty}dx~x\left[\frac{xe^x + e^x - 1}{\left(e^x-1\right)^2}\right]
Now, there's probably some clever way to evaluate the integral there in terms of polylogarithms. I could probably do it with some thinking, but I think for now I won't. Wolfram seems to be able to do it, but I can't evaluate the definite integral online, so instead I just evaluated it numerical, from x = 0.000000000001 (or so) to 49, which gives me 4.9348... . Dividing this by 12 \pi^2 gives a number pretty close to 1/24.
Hence, there's the second term of the series expansion. You can keep doing this game to generate higher order terms, though each one is probably going to give you more and more annoying-to-evaluate integrals.
So, to second order,
\ln Z_\pi \simeq (\nu \Omega \beta)T^4 \left( \frac{\pi^2}{90} - \frac{m_\pi^2}{24T^2} + \cdots~\right)
