I've been trying to solve the following definite integral $$ \int_0^\infty dx\, x^4\, \frac{e^{x^2+a}}{\left(e^{x^2+a}-1\right)^2}\quad , \qquad a>0\, . $$ However, so far I could not find any solution. I'd prefer an exact one but a good approximation would be also nice.
I know that the following relations hold: $$ \int_0^\infty x^d \frac{e^x}{\left(e^{x}-1\right)^2} = d\Gamma(d)\zeta(d)\\ \int_0^\infty x^{d-1} \frac{1}{e^{x}-1} = \Gamma(d)\zeta(d)\\ \int_0^\infty x^{d-1} \text{ln}\{1-e^{-x}\} = -\Gamma(d)\zeta(d+1)\, . $$ There are several ways I tried to solve this. Most of them are not worth mentioning. The most promising one is the following: I substituted $u\equiv x^2$ to get $$ \frac{1}{2}\int_0^\infty du\, u^{\frac{3}{2}}\, \frac{e^{u+a}}{\left(e^{u+a}-1\right)^2}\, . $$ Then, one can see that $$ \frac{1}{2}\int_0^\infty du\, u^{\frac{3}{2}}\, \frac{e^{u+a}}{\left(e^{u+a}-1\right)^2} = -\frac{1}{2}\frac{\partial}{\partial a}\int_0^\infty du\, u^{\frac{3}{2}}\, \frac{1}{e^{u+a}-1} $$ If one now substitutes $u^\prime\equiv u+a$ then the integration boundaries change from $1$ to $\infty$ such that one cannot use the above relations. For this reason, I considered an approximation for small $a$ around 0 by writing $$ \frac{\partial}{\partial a}\int_0^\infty du\, u^{\frac{3}{2}}\, \frac{1}{e^{u}(1+a)-1}\, . $$ Now, if one expands the integrand in $a$ up to infinite order, one gets $$ \frac{\partial}{\partial a}\int_0^\infty du\, u^{\frac{3}{2}}\, \frac{1}{e^{u}-1}\sum_{n=0}^\infty\left(-a\frac{e^u}{e^u-1}\right)^n\, . $$ I know that the integral of the expression in the sum is a Hypergeometric function so maybe one could use partial integration. It did not take me anywhere at least.
Does anyone happen to know the result of this or how I could solve it. I am also happy with a nice approximation.
Thanks in advance!
Denoting $t=e^{-a}$, substituting $x=y^{1/2}$ and using $z(1-z)^{-2}=\sum_{n=1}^\infty nz^n$, you get $\operatorname{Li}_{3/2}$: $$\ldots=\frac12\int_0^\infty\frac{y^{3/2}te^{-y}}{(1-te^{-y})^2}\,dy=\frac12\sum_{n=1}^\infty nt^n\int_0^\infty y^{3/2}e^{-ny}\,dy=\frac{3\sqrt\pi}{8}\sum_{n=1}^\infty\frac{t^n}{n^{3/2}}.$$