Evaluating a definite integral by introducing change of variables In a Physics class, I encountered the following integral: $$\int_{0}^{\infty} \left(\frac{x^2}{\exp\left(x-a\right)+1}-\frac{x^2}{\exp\left(x+a\right)+1}\right)\,\mathrm{d}x\text{,}$$
where $a$ is a constant. It is stated to yield: $$\int_{0}^{\infty}\left(\frac{x^2}{\exp\left(x-a\right)+1}-\frac{x^2}{\exp\left(x+a\right)+1}\right)\,\mathrm{d}x=\frac{1}{3}\left[\pi^2a+a^3\right]$$
I first introduced a change of variables: $$y:=x-a$$ $$z:=x^\prime+a$$ This gives: $$\int_{0}^{\infty}\frac{x^2}{\exp\left(x-a\right)+1}\,\mathrm{d}x-\int_{0}^{\infty}\frac{\left(x^\prime\right)^2}{\exp\left(x^\prime+a\right)+1}\,\mathrm{d}x^\prime=\int_{-a}^{\infty}\frac{\left(y+a\right)^2}{e^y + 1}\,\mathrm{d}y-\int_{a}^{\infty}\frac{\left(z-a\right)^2}{e^z+1}\,\mathrm{d}z$$ However, when introduced to a computer, the definite integrals are stated to have no solution. I actually expected to get integrals of the form: $$\int_{0}^{\infty}\frac{y}{e^y+1}\,\mathrm{d}y=\frac{\pi^2}{12}$$ What am I missing here?
 A: You are on the right track. Continue by renaming $z$ to $y$ and writing
\begin{align}
\int \limits_0^\infty x^2 \left(\frac{1}{\mathrm{e}^{x-a}+1} - \frac{1}{\mathrm{e}^{x+a}+1}\right) \, \mathrm{d} x &= \int \limits_{-a}^\infty \frac{(y+a)^2}{\mathrm{e}^y + 1} \, \mathrm{d} y - \int \limits_a^\infty \frac{(z-a)^2}{\mathrm{e}^z + 1} \, \mathrm{d} z \\
&= \int \limits_{-a}^a \frac{(y+a)^2}{\mathrm{e}^y + 1} \, \mathrm{d} y + \int \limits_a^\infty \frac{(y+a)^2 - (y-a)^2}{\mathrm{e}^{y} + 1} \, \mathrm{d} y \\
&= \int \limits_{-a}^0 \frac{(y+a)^2}{\mathrm{e}^y + 1} \, \mathrm{d} y + \int \limits_0^a \frac{(y+a)^2 - 4 a y}{\mathrm{e}^y + 1} \, \mathrm{d} y + \int \limits_0^\infty \frac{4ay}{\mathrm{e}^y + 1} \, \mathrm{d} y \\
&= \int \limits_{-a}^0 \frac{(a + y)^2}{\mathrm{e}^y + 1} \, \mathrm{d} y + \int \limits_0^a \frac{(a - y)^2}{\mathrm{e}^y + 1} \, \mathrm{d} y + \int \limits_0^\infty \frac{4ay}{\mathrm{e}^y + 1} \, \mathrm{d} y \, .
\end{align}
Now substitute $y \to - y$ in the first integral and add it to the second integral. Since you know the value of the third integral, this yields
\begin{align}
\int \limits_0^\infty x^2 \left(\frac{1}{\mathrm{e}^{x-a}+1} - \frac{1}{\mathrm{e}^{x+a}+1}\right) \, \mathrm{d} x &= \int \limits_0^a (a - y)^2 \left[\frac{1}{\mathrm{e}^{-y} + 1} + \frac{1}{\mathrm{e}^y + 1}\right] \, \mathrm{d} y + 4 a \int \limits_0^\infty \frac{y}{\mathrm{e}^y + 1} \, \mathrm{d} y \\
&= \int \limits_0^a (a-y)^2 \, \mathrm{d} y + 4 a \frac{\pi^2}{12} \\
&= \frac{a}{3}(\pi^2 + a^2)
\end{align}
for $a \in \mathbb{R}$ as desired. Note that we had to combine the two integrals over $y$ and $z$ in a particular way to obtain this result, which is why some computer programs may be unable to evaluate them individually.
A: If you encountered this integral in a statistical mechanics class, perhaps they wanted you to use the integral of the Fermi-Dirac distribution and properties of the polylogarithm function.
Since $$-\operatorname{Li}_{3}(-x) = \frac{1}{2} \int_{0}^{\infty} \frac{t^{2}}{e^{t}/x+1} \, \mathrm dt, \, \quad x \ge -1,$$ and the inversion formula for the trilogarithm function states that $$\operatorname{Li}_{3}(-x) - \operatorname{Li}_{3} \left(- \frac{1}{x} \right) = - \frac{1}{6} \, \ln^{3}(x)  - \frac{\pi^{2}}{6} \ln x \, , \quad x \notin (-1,0], \tag{1}$$
we have $$ \begin{align} \int_{0}^{\infty} \left(\frac{t^2}{e^{t-a}+1}-\frac{t^2}{e^{t+a}+1}\right)\,\mathrm dt &= - 2\operatorname{Li}_{3}(-e^{a})+ 2 \operatorname{Li}_{3} \left(-\frac{1}{e^{a}}\right) \\ &= 2 \left( \frac{1}{6} \, \ln^{3}(e^{a}) + \frac{\pi^{2}}{6} \ln(e^{a}) \right) \\ &= 2 \left(\frac{a^{3}}{6}+ \frac{ \pi^{2}a }{6} \right) \\ &= \frac{a}{3} \left(a^{2} + \pi^{2} \right). \end{align}$$

$(1)$ The inversion formula for the trilogarithm can be derived from the inversion formula of the dilogarithm by dividing both sides of the equation by $x$ and integrating.
