Evaluating $\int_0^\infty \sqrt{\frac{x}{e^x-1}}dx$ in terms of special functions

A partial answer.

I'm sure you already know the Lerch transcendent, a special function which may initially be defined as $$\Phi(z,s,a):=\sum_{k=0}^\infty\frac{z^k}{(a+k)^s}, \quad a>0,\Re s>1,|z|<1.$$ It admits the following integral representation, which you obtain by expanding the integrand as a powers series of $z$: $$ \Phi(z,s,a)=\int_0^{\infty}\frac{x^{s-1}e^{-ax}}{1-ze^{-x}}{\rm d}x. $$ By differentiation with respect to $z$, you get $$ \partial_z^r\Phi(z,s,a)=(-1)^r\int_0^{\infty}\frac{x^{s-1}e^{-(a+r)x}}{(1-ze^{-x})^{r+1}}{\rm d}x.\tag1 $$ This shows the level of complexity of your family of integrals: fractional calculus.

For example, from $(1)$, you have

$$ \int_0^\infty \sqrt\frac{x}{e^x-1}dx= \color{blue}{\left(-\partial_z\right)^{\large \! -1/2}\Phi\left(1,3/2,1\right). } $$

I would be very surprised to learn that we have more than the above expression.

Starting from the general $~\displaystyle\int_0^\infty\frac{x^a}{e^x-u}dx~=~\frac{\Gamma(a+1)\cdot\text{Li}_{a+1}(u)}u,~$ we can then

deduce, by way of repeated differentiation under the integral sign with regard to u,

that $~\displaystyle\int_0^\infty\frac{x^a}{(e^x-u)^b}dx~=~(-1)^b\cdot\Gamma(a+1)\cdot{\large\partial}^b_u~\bigg[\dfrac{\text{Li}_{a+1}(u)}u\bigg],~$ for all real $a>-1$

and $u\in\mathbb C$, with $b\in\mathbb N.~$ However, in this case, $~b=\dfrac12~$ is fractional, so, if you are

looking for a closed form, you will have to appeal to the concept of fractional calculus,

which has already been mentioned above, though in relation to a lesser-known special

function.