I would like to prove that for every real $\alpha>1$ we have $$\int_0^\infty \frac{x^\alpha}{\exp(x) -1} \, dx=\zeta(\alpha+1)\Gamma(\alpha+1).$$ Proof: Let $0<a<b$; then we have $$\int_a^b \frac{x^\alpha}{\exp(x) -1} \, dx=\int_a^b x^\alpha \exp(-x)\sum_{k=0}^\infty\exp(-kx)dx$$ as $1/(1-\exp(-x))=\sum_{k=0}^\infty\exp(-kx)$: this series converges uniformly on $[a,b]$ by the Weierstrass M-Test: $\exp(-kx)\leq\exp(-ka)$ on $[a,b]$ and $\exp(-a)<1$ since $a>0$, use geometric series. Thus we may exchange integral and series, giving $$\int_a^b \frac{x^\alpha}{\exp(x) -1} \, dx=\sum_{k=0}^\infty\int_a^b x^\alpha \exp(-(1+k)x) \, dx=\sum_{k=0}^\infty\frac{1}{(1+n)^\alpha}\int_{(1+k)a}^{(1+k)b} s^\alpha\exp(-s) \, ds$$ (Here was a glaring error, thanks to D. Fischer.) Now taking $a\downarrow 0$ and $b\uparrow\infty$ we see that $$\int_0^\infty \frac{x^\alpha}{\exp(x) -1} \, dx$$ exists, since $$\int_0^\infty s^\alpha \exp(-s) \, ds$$ does and is $=\Gamma(\alpha+1)$.
My questions: is the proof correct and complete? I am particularly interested in the technical points, such as whether the proof of existence of the integral is correct and whether the argument for the exchange of the integral and series is correct (the integral is taken to be a Riemann integral, of course).
EDIT: The question is now of course to fill the following gap: why can we exchange $\lim_{a\downarrow 0}\lim_{b\uparrow\infty}$ with the series?
Since the substitution $s = (1+k)x$ transforms
$$\int_a^b x^\alpha e^{-(1+k)x}\,dx$$
into
$$\frac{1}{(1+k)^{\alpha+1}} \int_{(1+k)a}^{(1+k)b} s^\alpha e^{-s}\,ds,$$
we can't pull the integral out of the sum as a constant, so the argument from a compact subinterval of $(0,\infty)$ does at least not work too easily.
The simplest way would of course be to use the Lebesgue integral, where the interchange of summation and integration is justified by one of the most fundamental theorems of the theory.
Within the framework of the Riemann integral, the easiest ways to prove the result that I see are
For the substitution, we could for example use $x = \frac{t}{1-t}$ and would obtain
$$\begin{align} \int_0^\infty \frac{x^\alpha}{e^x-1}\,dx &= \int_0^1 \frac{t^\alpha}{(1-t)^{\alpha+2} (e^{t/(1-t)}-1)}\,dt\\ &= \int_0^1 \frac{t^\alpha}{(1-t)^{\alpha+2}}\sum_{k=1}^\infty \exp\left(-\frac{kt}{1-t}\right)\,dt. \end{align}$$
The integrand
$$f(t) = \frac{t^\alpha}{(1-t)^{\alpha+2}(e^{t/(1-t)}-1)}$$
extends to a continuous function on $[0,1]$ with $f(0) = f(1) = 0$, as do the partial sums
$$f_N(t) = \frac{t^\alpha}{(1-t)^{\alpha+2}}\sum_{k=1}^N \exp\left(-\frac{kt}{1-t}\right),$$
and since the sequence is monotonic, the convergence is uniform, so we can interchange integration and summation. The one can substitute back ($s = kt/(1-t)$) and obtain
$$\int_0^\infty \frac{x^\alpha}{e^x-1}\,dx = \sum_{k=1}^\infty \frac{1}{k^{\alpha+1}}\int_0^\infty s^\alpha e^{-s}\,ds = \Gamma(\alpha+1)\zeta(\alpha+1).$$
For the version of the monotone convergence theorem that suffices here:
Let $f \colon (0,\infty) \to [0,\infty)$ a continuous function such that $\limsup\limits_{x\searrow 0} f(x) < \infty$ and
$$\sup_{0 < a < b < \infty} \int_a^b f(x)\,dx < \infty.$$
Let $(f_n)$ be a sequence of continuous functions $(0,\infty) \to [0,\infty)$ with $f_n(x) \leqslant f_{n+1}(x)$ and $f_n(x) \to f(x)$ for all $x\in (0,\infty)$. Then
$$\int_0^\infty f(x)\,dx;\qquad \int_0^\infty f_n(x)\,dx$$
exist as improper Riemann integrals, and
$$\int_0^\infty f(x)\,dx = \lim_{n\to\infty} \int_0^\infty f_n(x)\,dx.$$
For the proof, use that for every $\varepsilon > 0$ there exist $0 < a_\varepsilon < b_\varepsilon < \infty$ such that
$$\int_0^{a_\varepsilon} f(x)\,dx < \frac{\varepsilon}{3},\qquad \int_{b_\varepsilon}^\infty f(x)\,dx < \frac{\varepsilon}{3},$$
and that the convergence of the $f_n$ to $f$ is uniform on $[a_\varepsilon,b_\varepsilon]$, so for all large enough $n$
$$\int_{a_\varepsilon}^{b_\varepsilon} \lvert f(x)- f_n(x)\rvert\,dx < \frac{\varepsilon}{3}.$$
Then the interchange of summation and integration in
$$\int_0^\infty x^\alpha \sum_{k=1}^\infty e^{-kx}\,dx = \sum_{k=1}^\infty \int_0^\infty x^\alpha e^{-kx}\,dx$$
is justified.
The monotone convergence theorem for series is of course the special case of the monotone convergence theorem for Lebesgue integrals where the measure is the counting measure on $\mathbb{N}$, but that special case can of course be elementarily proved without mentioning $\sigma$-algebras, measures, or measurability.
Let $(a_n^{(k)})_{n\in\mathbb{N}}$ a sequence of non-negative real numbers for each $k\in \mathbb{N}$, and suppose for all $n,k$ we have $a_n^{(k)} \leqslant a_n^{(k+1)}$. Let $a_n = \lim\limits_{k\to\infty} a_n^{(k)}$. Then
$$\sum_{n=0}^\infty a_n = \lim_{k\to\infty} \sum_{n=0}^\infty a_n^{(k)}.$$
Evidently $a_n^{(k)}\leqslant a_n$ for all $n,k$, whence
$$S_k := \sum_{n=0}^\infty a_n^{(k)} \leqslant \sum_{n=0}^\infty a_n =: S$$
for all $k$, and $S_k \leqslant S_{k+1}$. Thus $\lim\limits_{k\to\infty} S_k \leqslant S$. On the other hand, let $K < S$. Then there is an $N\in\mathbb{N}$ with
$$\sum_{n=0}^N a_n > K,$$
and we can find a $k\in\mathbb{N}$ such that
$$a_n^{(k)} > a_n - \frac{1}{2(N+1)}\left(\sum_{n=0}^N a_n - K\right)$$
for $0 \leqslant n \leqslant N$ by monotonicity. Then
$$S_k \geqslant \sum_{n=0}^N a_n^{(k)} > \sum_{n=0}^N\left(a_n - \frac{1}{2(N+1)}\left(\sum_{n=0}^N a_n - K\right)\right) = \frac{1}{2} \left(K + \sum_{n=0}^N a_n\right) > K,$$
whence $\lim\limits_{k\to\infty} S_k > K$. Since that holds for all $K < S$, the proposition $\lim\limits_{k\to\infty} S_k = S$ follows.
Then we can choose monotonic sequences $a_m \searrow 0$, and $b_m \nearrow\infty$, and applying the monotone convergence theorem to
$$S_m := \int_{a_m}^{b_m} \frac{x^\alpha}{\exp(x) -1} \, dx=\sum_{k=0}^\infty\int_{a_m}^{b_m} x^\alpha \exp(-(1+k)x) \, dx=\sum_{k=0}^\infty\frac{1}{(1+k)^\alpha}\int_{(1+k)a_m}^{(1+k)b_m} s^\alpha\exp(-s) \, ds$$
with
$$\lim_{m\to\infty} \frac{1}{(1+k)^\alpha}\int_{(1+k)a_m}^{(1+k)b_m} s^\alpha\exp(-s) \, ds = \frac{\Gamma(\alpha+1)}{(1+k)^{\alpha+1}}$$
then yields the result.