When describing black body radiation in Physics, one comes across an integral of the form
$$I(a)=\int_0^\infty\frac{x^3}{e^{ax}-1}dx \tag1$$
for $a>0$. In our lecture, we were given that this integral evaluates to
$$I(a)=\frac{\pi^4}{15a^4}$$
without a derivation of this result.
While reviewing my course notes over the Christmas break, I came across this integral and have tried to evaluate it myself, but I have not been successful so far.
My first step, of course, was to substitute $u:=ax$, which gives us
$$I(a) = \frac{1}{a^4}\int_0^\infty\frac{u^3}{e^u-1}du \tag2$$
but beyond that, I'm stuck.
I've tried differentiation under the integral sign on the original integral in $(1)$, which lead to the integral
$$\int_0^\infty\frac{\log^4 u}{(u-1)^2}du$$
which I also no have no clue how to tackle.
I've also thought about using contour integration, but I'm unsure as to what would be an appropriate contour to consider. I know that the integrand has singularities at $z=2\pi ik$ for $k\in\mathbb{Z},k\neq 0$ but I'm not even sure how one would go about calculating the residues at these singularities due to the peculiar form of the denominator, which is unlike any I have previously encountered.
How would one go about deriving that
$$\int_0^\infty\frac{x^3}{e^{ax}-1}dx = \frac{\pi^4}{15a^4}$$
is true?