Evaluation of Bose-Einstein and Fermi-Dirac Integrals Can anyone help me with the following integral, which has arisen from my study of cosmology?
$$\int_0^\infty\frac{p^2dp}{e^{(p\pm\mu)/T}+ 1}$$
I have tried substitution and using standard integral tables in Gradshteyn and Ryzhik, but the change of variables alters the lower limit on the integral, so I'm not sure if standard definite integrals can then be applied.
I have also tried integrating by parts, but this just leads me to further similar integrals that I cannot evaluate.
I would basically like to know if the integral can be evaluated exactly, or if it is necessary to use a series expansion for the denominator.
Update
I have found an exact relation involving the original integral: 
$$\int_0^\infty\left(\frac{1}{e^{(p-\mu)/T}+1}-\frac{1}{e^{(p+\mu)/T}+1}\right)p^2\;dp=\frac{T^3}{3}\left[\pi^2\left(\frac{\mu}{T}\right)+\left(\frac{\mu}{T}\right)^3\right]$$  
However, I am not sure how this was evaluated either, but maybe it can help?
 A: This is a derivation of the integral result provided by the OP in the Update.  The integral is
$$\begin{align}\int_0^{\infty} dp \, p^2 \left (\frac1{e^{(p-\mu)/T}+1} - \frac1{e^{(p+\mu)/T}+1} \right ) &= 2 \sinh{\frac{\mu}{T}} \int_0^{\infty} dp \;\frac{p^2 \, e^{p/T}}{e^{2 p/T}+ 2 \cosh{\frac{\mu}{T}}\; e^{p/T} + 1} \\ &=  2 T^3 \sinh{\frac{\mu}{T}} \int_1^{\infty} dx \; \frac{\log^2{x}}{x^2+2\cosh{\frac{\mu}{T}}\; x+1}\\ &= T^3 \sinh{\frac{\mu}{T}} \int_0^{\infty} dx \; \frac{\log^2{x}}{x^2+2\cosh{\frac{\mu}{T}}\; x+1} \end{align}$$
In the second line, I subbed $p=T \log{x}$.
The last integral may be evaluated via contour integration by considering the following contour integral:
$$\oint_C dz \frac{\log^3{z}}{z^2+2\cosh{\frac{\mu}{T}} \; z+1} $$
where $C$ is a keyhole contour of outer radius $R$ and inner radius $\epsilon$ about the positive real axis, i.e., so that $\arg{z} \in [0,2 \pi]$.  The contour integral is, as $R \to \infty$ and $\epsilon \to 0$,
$$\int_0^{\infty} dx \frac{\log^3{x} - (\log{x}+i 2 \pi)^3}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} $$
which is equal to
$$-i 6 \pi \int_0^{\infty} dx \frac{\log^2{x}}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} + 12 \pi^2 \int_0^{\infty} dx \frac{\log{x}}{x^2+2\cosh{\frac{\mu}{T}} \; x+1}\\+ i 8 \pi^3 \int_0^{\infty}  \frac{dx}{x^2+2\cosh{\frac{\mu}{T}} \; x+1}$$
Note that the second integral vanishes.  You can show this by subbing $x = 1/y$ and seeing that the integral is equal to the negative of itself.  
On the other hand, the contour integral is also equal to the sum of the residues at the poles $z_+=-e^{\mu/T}$ and $z_-=-e^{-\mu/T}$.  The negative signs will be treated as factors of $e^{i \pi}$ because of the branch of the log we are using as defined by the contour $C$.  Thus, by the residue theorem,
$$-i 6 \pi \int_0^{\infty} dx \frac{\log^2{x}}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} + i 8 \pi^3 \int_0^{\infty}  \frac{dx}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} = i 2 \pi \left [\frac{\left (\frac{\mu}{T} + i \pi \right )^3}{-2 \sinh{\frac{\mu}{T}}} + \frac{\left (-\frac{\mu}{T} + i \pi \right )^3}{2 \sinh{\frac{\mu}{T}}} \right ]$$
We're almost done...except that we need the second integral on the left.  The good news is that we can evaluate it in precisely the same manner as we did the original integral, but with a factor of $\log{z}$ rather than $\log^3{z}$ in the numerator of the integrand.  Thus,
$$-i 2 \pi \int_0^{\infty}  \frac{dx}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} = i 2 \pi \left [\frac{\frac{\mu}{T} + i \pi }{-2 \sinh{\frac{\mu}{T}}} + \frac{-\frac{\mu}{T} + i \pi }{2 \sinh{\frac{\mu}{T}}} \right ]$$
or
$$\int_0^{\infty}  \frac{dx}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} = \frac{\mu}{T \sinh{\frac{\mu}{T}}} $$
Thus, we have all the pieces in place to evaluate the integral.  We find that, upon expansion of the RHS and combining like terms, we get
$$\int_0^{\infty} dx \frac{\log^2{x}}{x^2+2\cosh{\frac{\mu}{T}} \; x+1} = \frac1{3 \sinh{\frac{\mu}{T}}} \left [\left (\frac{\mu}{T} \right )^3 + \pi^2 \frac{\mu}{T} \right ]$$
Thus, our original integral is

$$\int_0^{\infty} dp \, p^2 \left (\frac1{e^{(p-\mu)/T}+1} - \frac1{e^{(p+\mu)/T}+1} \right ) = \frac{T^3}{3} \left [\left (\frac{\mu}{T} \right )^3 + \pi^2 \frac{\mu}{T} \right ]$$

as asserted.
A: Just for fun (After all, integration is just art).
Suppose you set $(p - \mu)/T = z$ so d$z = \text{d}p/T$ and $p = zT + \mu$.
In this way you get:
$$\int_{-\mu/T}^{\infty} T\frac{(zT + \mu)^2}{e^z \pm 1}\ \text{d}z$$
you may collect $e^z$, getting
$$T\int_{-\mu/T}^{\infty}\frac{(zT + \mu)^2}{e^z (1 \pm e^{-z})}\ \text{d}z$$
Now usually when we deal with integral of this kind, we pull out a geometric series from the term
$$\frac{1}{1 \pm e^{-z}}$$
since $z$ is usually in the range $[0, +\infty)$. In this case the range is different, starting from a negative $\mu/T$ which could forbid us to generate the geometric series. We may have a try the same. The trick is the same either for $+$ and for $-$ so I'll give you the $-$ treatise.
$$\frac{1}{1 - e^{-z}} = \sum_{k = 0}^{+\infty} (e^{-z})^k$$
thence
$$T\int_{-\mu/T}^{\infty} (zT + \mu)^2 e^{-z}\sum_{k = 0}^{+\infty} (e^{-zk})\ \text{d}z = T\sum_{k = 0}^{+\infty}\int_{-\mu/T}^{\infty} (zT + \mu)^2 e^{-z}(e^{-zk})\ \text{d}z$$
I think there is no close form for such an integration, but you can always try to split the square and see what happens..
$$T\sum_{k = 0}^{+\infty}\int_{-\mu/T}^{\infty} \left(z^2T^2 e^{-z(1+k)} + \mu^2 e^{-z(1+k)} + 2zT\mu e^{-z(1+k)}\right)\ \text{d}z$$
Remark
Due to the lower extrema of the integral, I'm not sure about the Geometric Series trick. I advice you to search for "Fermi-Dirac" integrals or "Bose-Einstein" integrals, since it appears exactly like one of them.
Let's proceed
Splitting the integrals into the three terms (leaving apart the series for the moment) we have the three integrals (valid for $\Re(k) > 0$) as follows:
$$\int_{-\mu/T}^{+\infty} z^2T^2e^{-z(1+k)}\ \text{d}z = T^2\frac{e^{\frac{\mu}{T}(1+k)}\left(2 - \frac{\mu}{T}(1+k)\left(2 - \frac{\mu}{T}- \frac{\mu}{T}k\right)\right)}{(1+k)^3}$$
$$\int_{-\mu/T}^{+\infty} \mu^2 e^{-z(1+k)}\ \text{d}z = \mu^2 \frac{e^{\frac{\mu}{T}(1+k)}}{1+k}$$
$$\int_{-\mu/T}^{+\infty} 2z\mu Te^{-z(1+k)}\ \text{d}z = 2T\mu\frac{e^{\frac{\mu}{T}(1+k)} \left(1 - \frac{\mu}{T} -\frac{\mu}{T}k\right)}{(1+k)^2}$$
For each of these integrals, you have a series. With the help of the PoliLog special functions, you could evaluate those three series:
$$T\sum_{k = 0}^{+\infty}T^2\frac{e^{\frac{\mu}{T}(1+k)}\left(2 - \frac{\mu}{T}(1+k)\left(2 - \frac{\mu}{T}- \frac{\mu}{T}k\right)\right)}{(1+k)^3}$$
$$T\sum_{k = 0}^{+\infty} \mu^2 \frac{e^{\frac{\mu}{T}(1+k)}}{1+k}$$
$$T\sum_{k = 0}^{+\infty} 2T\mu\frac{e^{\frac{\mu}{T}(1+k)} \left(1 - \frac{\mu}{T} -\frac{\mu}{T}k\right)}{(1+k)^2}$$
Summation of Series - Mathematica will help
I'm going now to write the summation of the three series saw below.
$$T\sum_{k = 0}^{+\infty}T^2\frac{e^{\frac{\mu}{T}(1+k)}\left(2 - \frac{\mu}{T}(1+k)\left(2 - \frac{\mu}{T}- \frac{\mu}{T}k\right)\right)}{(1+k)^3} = \frac{T^3}{8}\left(\frac{\mu^2}{T^2}e^{2\mu /T}\ _5F_4\left[\{ 2, 2, 2, 2, 2\}, \{1, 3, 3, 3\}e^{\mu /T}\right] - 16\text{PolyLog}[2, e^{\mu /T}] + 16\frac{\mu^2}{T^2}\text{PolyLog}[2, e^{\mu /T}] + 16\text{PolyLog}[3, e^{\mu /T}] - 8\frac{\mu^2}{T^2}\text{PolyLog}[3, e^{\mu /T}]\right)$$
$$T\sum_{k = 0}^{+\infty} \mu^2 \frac{e^{\frac{\mu}{T}(1+k)}}{1+k} = T\mu^2\ \ln\left(1 - e^{\mu /T}\right)$$
$$T\sum_{k = 0}^{+\infty} 2T\mu\frac{e^{\frac{\mu}{T}(1+k)} \left(1 - \frac{\mu}{T} -\frac{\mu}{T}k\right)}{(1+k)^2} = 2T^2\mu\left(\frac{\mu}{T}\ln\left(1- e^{\mu}{T}\right) + \text{PolyLog}[2, e^{\mu /T}]\right)$$
End of the Hell
Thence the solution by series of your integral becomes:
$$
\frac{T^3}{8}\left(\frac{\mu^2}{T^2}e^{2\mu /T}\ _5F_4\left[\{ 2, 2, 2, 2, 2\}, \{1, 3, 3, 3\}e^{\mu /T}\right] - 16\text{PolyLog}[2, e^{\mu /T}] + 16\frac{\mu^2}{T^2}\text{PolyLog}[2, e^{\mu /T}] + 16\text{PolyLog}[3, e^{\mu /T}] - 8\frac{\mu^2}{T^2}\text{PolyLog}[3, e^{\mu /T}]\right) + T\mu^2\ \ln\left(1 - e^{\mu /T}\right) + 2T^2\mu\left(\frac{\mu}{T}\ln\left(1- e^{\mu}{T}\right) + \text{PolyLog}[2, e^{\mu /T}]\right)
$$
More on special function involved
Here for The generalized Hypergeometric Function $_qF_p$
https://reference.wolfram.com/language/ref/HypergeometricPFQ.html
Here for the PolyLogarithm Function
https://en.wikipedia.org/wiki/Polylogarithm
http://mathworld.wolfram.com/Polylogarithm.html
A: By the change of variable $\dfrac{p-\mu}T=t$, the integral is
$$T\int_{-\mu/T}^\infty \frac{(Tt+\mu)^2dt}{e^{t/T}\pm1}.$$ 
Because of the variable lower bound, you can consider it as an indefinite integral.
Then Wolfram gives no good news: the terms in $t^2$ and $t$ do imply expressions with $Li_2$ and $Li_3$, and no possible cancellation.
