Closed form solution to normal fluid density integral in the two fluid model 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\, .
$$
This is the density of the normal part of a superfluid. 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}\, ,
$$
where the fraction part is the Bose distribution.
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.
Edit: In fact the above integral is already an approximation. The original integral was
$$
\int_0^\infty dx\, x^4\, \frac{e^{\sqrt{x^4+2 x^2}/Tp}}{\left(e^{\sqrt{x^4+2 x^2}/Tp}-1\right)^2}\, ,\quad Tp = \frac{T}{Un}
$$
I'd be happy if anyone could point me to a solution if there is one...
 A: Let define $f(a)=\int_0^{\infty}\frac{x^4e^{x^2+a}}{(e^{x^2+a}-1)^2}$, then $f(a)=-\frac{\partial}{\partial a}\int_0^{\infty}\frac{x^4}{e^{x^2+a}-1}$. Last integral is calculated exactly as $\frac {3}{8}\sqrt{\pi}Li_{\frac{5}{2}}(e^{-a})$, so $f(a)=\frac {3}{8}\sqrt{\pi}Li_{\frac{3}{2}}(e^{-a})$, where $Li_n(z)$ is the polylogarithm function. In Figure 1 shown function $f(a)$ (solid line)  with a numerical calculation (points)

A: Where do you get your expression for the normal fluid density? For superfluid helium the usual expression for  is $\rho_n$ at low superfluid velicity is
$$
\rho_n|_{v_s\to 0}= \frac 1{3T} \int \frac{d^3p}{(2\pi)^3} \frac {p^2 e^{E(p)/T}}{(e^{E(p)/T}-1)^2}
$$
Here, for $|p|$ below the roton minimum, we can take  $E(p) = c_{\rm sound} |p|$. The integral looks like yours, but I don't see where your  $E(p)$ comes from. Are you just considering the rotons to  get $E(x)= x^2+a$?
There is a detailed discussion of both the phonon and the roton  roton contribution  here.
