[proofplan]
We start from the integral representation of the Gamma function, $\Gamma(s) = \int_0^\infty t^{s-1} e^{-t}\, d\mathcal{L}^1(t)$, and substitute $t \mapsto nt$ to express $n^{-s}\Gamma(s)$ as an integral. Summing over $n \geq 1$ and interchanging the sum and integral (justified by Tonelli's theorem for $\operatorname{Re}(s) > 1$) produces the geometric series $\sum_{n=1}^\infty e^{-nt} = (e^t - 1)^{-1}$, yielding the desired formula.
[/proofplan]
[step:Express $n^{-s}\Gamma(s)$ as an integral via substitution in the Gamma function]
The Gamma function for $\operatorname{Re}(s) > 1$ is defined by
\begin{align*}
\Gamma(s) &= \int_0^\infty t^{s-1} e^{-t}\, d\mathcal{L}^1(t).
\end{align*}
For each integer $n \geq 1$, perform the substitution $t \mapsto nt$ (i.e., set $u = nt$, so $t = u/n$ and $d\mathcal{L}^1(t) = n^{-1}\, d\mathcal{L}^1(u)$, and the domain $t \in (0,\infty)$ maps to $u \in (0,\infty)$). Then
\begin{align*}
\Gamma(s) &= \int_0^\infty \left(\frac{u}{n}\right)^{s-1} e^{-u}\, \frac{d\mathcal{L}^1(u)}{n} = n^{-s} \int_0^\infty u^{s-1} e^{-u}\, d\mathcal{L}^1(u).
\end{align*}
Wait — this is circular. Instead, substitute directly in the original integral: for $n \geq 1$, substitute $u = nt$ in the integral $\int_0^\infty u^{s-1} e^{-u}\, d\mathcal{L}^1(u) = \Gamma(s)$. Setting $u = nt$, $d\mathcal{L}^1(u) = n\, d\mathcal{L}^1(t)$, the domain remains $(0, \infty)$:
\begin{align*}
\Gamma(s) &= \int_0^\infty (nt)^{s-1} e^{-nt}\, n\, d\mathcal{L}^1(t) = n^s \int_0^\infty t^{s-1} e^{-nt}\, d\mathcal{L}^1(t).
\end{align*}
Dividing both sides by $n^s$:
\begin{align*}
\frac{\Gamma(s)}{n^s} &= \int_0^\infty t^{s-1} e^{-nt}\, d\mathcal{L}^1(t).
\end{align*}
[/step]
[step:Sum over $n$ and interchange sum and integral via Tonelli's theorem]
Summing over $n = 1, 2, 3, \ldots$:
\begin{align*}
\Gamma(s) \sum_{n=1}^\infty \frac{1}{n^s} &= \sum_{n=1}^\infty \int_0^\infty t^{s-1} e^{-nt}\, d\mathcal{L}^1(t).
\end{align*}
The left-hand side is $\Gamma(s)\,\zeta(s)$ for $\operatorname{Re}(s) > 1$ (where the Dirichlet series converges absolutely). To interchange the sum and integral on the right-hand side, we apply Tonelli's theorem. We must verify that the summands are non-negative measurable functions or that the double integral of the absolute value is finite.
Write $s = \sigma + i\tau$ with $\sigma = \operatorname{Re}(s) > 1$. The absolute value of the integrand is
\begin{align*}
|t^{s-1} e^{-nt}| &= t^{\sigma - 1} e^{-nt}
\end{align*}
for $t > 0$. This is non-negative, and
\begin{align*}
\sum_{n=1}^\infty \int_0^\infty t^{\sigma-1} e^{-nt}\, d\mathcal{L}^1(t) &= \sum_{n=1}^\infty \frac{\Gamma(\sigma)}{n^\sigma} = \Gamma(\sigma)\,\zeta(\sigma) < \infty
\end{align*}
since $\sigma > 1$. The finiteness of this sum of integrals (with non-negative integrands) verifies the hypothesis of Tonelli's theorem on the product measure space $(\mathbb{N}, \text{counting}) \times ((0,\infty), \mathcal{L}^1)$. Therefore we may interchange:
\begin{align*}
\sum_{n=1}^\infty \int_0^\infty t^{s-1} e^{-nt}\, d\mathcal{L}^1(t) &= \int_0^\infty t^{s-1} \sum_{n=1}^\infty e^{-nt}\, d\mathcal{L}^1(t).
\end{align*}
[guided]
Tonelli's theorem applies to non-negative measurable functions on a product of $\sigma$-finite measure spaces, allowing interchange of the order of integration (equivalently, sum and integral) without any integrability hypothesis beyond measurability and non-negativity. Here the measure spaces are $(\mathbb{N}, 2^{\mathbb{N}}, \text{counting measure})$ (which is $\sigma$-finite) and $((0,\infty), \mathcal{B}((0,\infty)), \mathcal{L}^1)$ (also $\sigma$-finite). The function $(n, t) \mapsto t^{\sigma-1} e^{-nt}$ is non-negative and measurable.
After applying Tonelli to the absolute values and confirming finiteness, we can also apply Fubini's theorem to the original complex-valued integrand (since $|t^{s-1}e^{-nt}| = t^{\sigma-1}e^{-nt}$ and the double integral of absolute values is finite). This justifies the interchange for the complex-valued function.
[/guided]
[/step]
[step:Evaluate the geometric series and conclude]
For $t > 0$, the geometric series converges:
\begin{align*}
\sum_{n=1}^\infty e^{-nt} &= \sum_{n=1}^\infty (e^{-t})^n = \frac{e^{-t}}{1 - e^{-t}} = \frac{1}{e^t - 1}.
\end{align*}
The ratio $e^{-t}$ satisfies $|e^{-t}| = e^{-t} < 1$ for $t > 0$, so the geometric series formula applies. Substituting into the integral:
\begin{align*}
\Gamma(s)\,\zeta(s) &= \int_0^\infty t^{s-1} \cdot \frac{1}{e^t - 1}\, d\mathcal{L}^1(t) = \int_0^\infty \frac{t^{s-1}}{e^t - 1}\, d\mathcal{L}^1(t).
\end{align*}
This is the desired integral representation, valid for $\operatorname{Re}(s) > 1$.
[guided]
Let us verify the integral on the right-hand side converges for $\operatorname{Re}(s) > 1$. Near $t = 0$, we have $e^t - 1 \sim t$ as $t \to 0^+$, so $t^{s-1}/(e^t - 1) \sim t^{s-2}$. The integral $\int_0^1 t^{\sigma - 2}\, d\mathcal{L}^1(t)$ converges iff $\sigma - 2 > -1$, i.e., $\sigma > 1$. Near $t = \infty$, the integrand decays like $t^{\sigma-1} e^{-t}$, which is integrable for any $\sigma$. So the integral converges precisely for $\operatorname{Re}(s) > 1$, matching the half-plane of absolute convergence of the Dirichlet series for $\zeta(s)$.
[/guided]
[/step]