[proofplan]
We follow the classical approach via the Jacobi theta function. Define $\theta(t) = \sum_{n=-\infty}^{\infty} e^{-\pi n^2 t}$ for $t > 0$, which satisfies the functional equation $\theta(1/t) = \sqrt{t}\,\theta(t)$. We express $\pi^{-s/2}\Gamma(s/2)\zeta(s)$ as a Mellin transform of $(\theta(t) - 1)/2$, split the integral at $t = 1$, and use the theta functional equation to transform the integral over $(0,1)$ into one over $(1,\infty)$. The resulting expression is manifestly meromorphic on all of $\mathbb{C}$ with a simple pole only at $s = 1$ (and $s = 0$, which is cancelled), and symmetric under $s \mapsto 1 - s$.
[/proofplan]
[step:Establish the integral representation of $\pi^{-s/2}\Gamma(s/2)\zeta(s)$ for $\operatorname{Re}(s) > 1$]
Define the completed zeta function
\begin{align*}
\xi(s) := \pi^{-s/2}\,\Gamma\!\left(\frac{s}{2}\right)\zeta(s).
\end{align*}
From the [Integral Representation of the Zeta Function](/theorems/3373), we have $\Gamma(s)\zeta(s) = \int_0^\infty t^{s-1}(e^t - 1)^{-1}\, d\mathcal{L}^1(t)$. We derive an alternative representation more suited to the theta function.
Starting from the Gamma function, for $\operatorname{Re}(s) > 0$:
\begin{align*}
\Gamma\!\left(\frac{s}{2}\right) &= \int_0^\infty t^{s/2 - 1} e^{-t}\, d\mathcal{L}^1(t).
\end{align*}
Substitute $t = \pi n^2 u$ for each integer $n \geq 1$ (so $d\mathcal{L}^1(t) = \pi n^2\, d\mathcal{L}^1(u)$, domain unchanged):
\begin{align*}
\Gamma\!\left(\frac{s}{2}\right) &= \int_0^\infty (\pi n^2 u)^{s/2 - 1} e^{-\pi n^2 u}\, \pi n^2\, d\mathcal{L}^1(u) = \pi^{s/2} n^s \int_0^\infty u^{s/2 - 1} e^{-\pi n^2 u}\, d\mathcal{L}^1(u).
\end{align*}
Dividing by $\pi^{s/2} n^s$:
\begin{align*}
\frac{\Gamma(s/2)}{\pi^{s/2} n^s} &= \int_0^\infty u^{s/2 - 1} e^{-\pi n^2 u}\, d\mathcal{L}^1(u).
\end{align*}
Summing over $n = 1, 2, 3, \ldots$ (interchange justified by Tonelli's theorem as in the proof of theorem 3373, since $\operatorname{Re}(s) > 1$ ensures absolute convergence of $\sum n^{-\sigma}$):
\begin{align*}
\pi^{-s/2}\,\Gamma\!\left(\frac{s}{2}\right)\zeta(s) &= \int_0^\infty u^{s/2 - 1} \sum_{n=1}^\infty e^{-\pi n^2 u}\, d\mathcal{L}^1(u) = \int_0^\infty u^{s/2 - 1}\,\psi(u)\, d\mathcal{L}^1(u),
\end{align*}
where we define the auxiliary function
\begin{align*}
\psi: (0, \infty) &\to \mathbb{R} \\
u &\mapsto \sum_{n=1}^\infty e^{-\pi n^2 u} = \frac{\theta(u) - 1}{2},
\end{align*}
and $\theta(u) := \sum_{n=-\infty}^{\infty} e^{-\pi n^2 u}$ is the Jacobi theta function.
[/step]
[step:Derive the functional equation of the Jacobi theta function]
The Jacobi theta function $\theta: (0,\infty) \to \mathbb{R}$ defined by $\theta(t) = \sum_{n=-\infty}^{\infty} e^{-\pi n^2 t}$ satisfies
\begin{align*}
\theta(1/t) &= \sqrt{t}\,\theta(t) \quad \text{for all } t > 0.
\end{align*}
This follows from the [Poisson Summation Formula](/theorems/???). For the Gaussian $g_t: \mathbb{R} \to \mathbb{R}$, $g_t(x) = e^{-\pi t x^2}$, its Fourier transform (under the symmetric normalisation) is
\begin{align*}
\hat{g}_t(\xi) &= \frac{1}{\sqrt{t}}\, e^{-\pi \xi^2 / t} = \frac{1}{\sqrt{t}}\, g_{1/t}(\xi).
\end{align*}
The Poisson summation formula states $\sum_{n \in \mathbb{Z}} g_t(n) = \sum_{n \in \mathbb{Z}} \hat{g}_t(n)$ (valid since $g_t \in \mathcal{S}(\mathbb{R})$), giving
\begin{align*}
\theta(t) = \sum_{n \in \mathbb{Z}} e^{-\pi n^2 t} = \sum_{n \in \mathbb{Z}} \frac{1}{\sqrt{t}} e^{-\pi n^2 / t} = \frac{1}{\sqrt{t}}\,\theta(1/t).
\end{align*}
Rearranging: $\theta(1/t) = \sqrt{t}\,\theta(t)$.
In terms of $\psi(t) = (\theta(t) - 1)/2$, this becomes
\begin{align*}
\psi(1/t) &= \frac{\theta(1/t) - 1}{2} = \frac{\sqrt{t}\,\theta(t) - 1}{2} = \frac{\sqrt{t}(2\psi(t) + 1) - 1}{2} = \sqrt{t}\,\psi(t) + \frac{\sqrt{t} - 1}{2}.
\end{align*}
[guided]
The Poisson summation formula is the engine driving the functional equation. For $f \in \mathcal{S}(\mathbb{R})$, it states $\sum_{n \in \mathbb{Z}} f(n) = \sum_{n \in \mathbb{Z}} \hat{f}(n)$. We need to check that $g_t(x) = e^{-\pi t x^2}$ is in the Schwartz class: it is smooth and decays (together with all derivatives) faster than any polynomial, so $g_t \in \mathcal{S}(\mathbb{R})$.
To compute $\hat{g}_t$: by definition (symmetric normalisation),
\begin{align*}
\hat{g}_t(\xi) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\pi t x^2} e^{-i\xi x}\, d\mathcal{L}^1(x).
\end{align*}
Actually, with the convention $\hat{f}(\xi) = \int_{\mathbb{R}} f(x) e^{-2\pi i \xi x}\, d\mathcal{L}^1(x)$ (which is more natural for Poisson summation), we get:
\begin{align*}
\hat{g}_t(\xi) = \int_{-\infty}^{\infty} e^{-\pi t x^2} e^{-2\pi i \xi x}\, d\mathcal{L}^1(x) = \frac{1}{\sqrt{t}} e^{-\pi \xi^2/t}
\end{align*}
by completing the square in the exponent: $-\pi t x^2 - 2\pi i \xi x = -\pi t(x + i\xi/t)^2 - \pi\xi^2/t$, and the Gaussian integral $\int e^{-\pi t w^2}\, d\mathcal{L}^1(w) = 1/\sqrt{t}$ (after contour shift, justified since the integrand is entire and decays in vertical strips).
With this convention of Fourier transform, Poisson summation reads $\sum_{n \in \mathbb{Z}} f(n) = \sum_{n \in \mathbb{Z}} \hat{f}(n)$, giving $\theta(t) = t^{-1/2}\theta(1/t)$, or equivalently $\theta(1/t) = \sqrt{t}\,\theta(t)$.
[/guided]
[/step]
[step:Split the integral at $t = 1$ and use the theta functional equation to obtain a meromorphic expression]
Starting from $\xi(s) = \int_0^\infty t^{s/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t)$ (valid for $\operatorname{Re}(s) > 1$), split at $t = 1$:
\begin{align*}
\xi(s) &= \int_0^1 t^{s/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t) + \int_1^\infty t^{s/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t).
\end{align*}
In the first integral, substitute $t \mapsto 1/t$ (so $d\mathcal{L}^1(t) \mapsto t^{-2}\, d\mathcal{L}^1(t)$, and the domain $(0,1)$ maps to $(1,\infty)$):
\begin{align*}
\int_0^1 t^{s/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t) &= \int_1^\infty t^{-(s/2 - 1)}\,\psi(1/t)\, t^{-2}\, d\mathcal{L}^1(t) = \int_1^\infty t^{-s/2 - 1}\,\psi(1/t)\, d\mathcal{L}^1(t).
\end{align*}
Using the functional equation $\psi(1/t) = \sqrt{t}\,\psi(t) + (\sqrt{t} - 1)/2$:
\begin{align*}
\int_1^\infty t^{-s/2 - 1}\left[\sqrt{t}\,\psi(t) + \frac{\sqrt{t} - 1}{2}\right] d\mathcal{L}^1(t) &= \int_1^\infty t^{(1-s)/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t) + \frac{1}{2}\int_1^\infty \left(t^{-s/2 - 1/2} - t^{-s/2 - 1}\right) d\mathcal{L}^1(t).
\end{align*}
The last integral evaluates explicitly: for $\operatorname{Re}(s) > 1$,
\begin{align*}
\int_1^\infty t^{-s/2 - 1/2}\, d\mathcal{L}^1(t) &= \frac{1}{s/2 - 1/2} = \frac{2}{s - 1}, \\
\int_1^\infty t^{-s/2 - 1}\, d\mathcal{L}^1(t) &= \frac{1}{s/2} = \frac{2}{s}.
\end{align*}
Therefore:
\begin{align*}
\xi(s) &= \int_1^\infty t^{(1-s)/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t) + \int_1^\infty t^{s/2 - 1}\,\psi(t)\, d\mathcal{L}^1(t) + \frac{1}{2}\left(\frac{2}{s-1} - \frac{2}{s}\right).
\end{align*}
Simplifying the rational term:
\begin{align*}
\frac{1}{s-1} - \frac{1}{s} &= \frac{1}{s(s-1)}.
\end{align*}
So the final formula is
\begin{align*}
\xi(s) &= \frac{1}{s(s-1)} + \int_1^\infty \left(t^{s/2 - 1} + t^{(1-s)/2 - 1}\right)\psi(t)\, d\mathcal{L}^1(t).
\end{align*}
[guided]
The substitution $t \mapsto 1/t$ is the standard Mellin-transform trick for exploiting the symmetry of the theta function. After substitution, the functional equation of $\psi$ introduces two parts: a "transformed" integral (with $s$ replaced by $1-s$) and explicit rational terms from the correction $(\sqrt{t}-1)/2$.
The computation of $\int_1^\infty t^{-\alpha}\, d\mathcal{L}^1(t)$ is elementary: for $\operatorname{Re}(\alpha) > 1$, $\int_1^\infty t^{-\alpha}\, d\mathcal{L}^1(t) = [t^{1-\alpha}/(1-\alpha)]_1^\infty = 1/(\alpha - 1)$. We apply this with $\alpha = s/2 + 1/2$ (giving $1/(s/2 - 1/2) = 2/(s-1)$) and $\alpha = s/2 + 1$ (giving $1/(s/2) = 2/s$).
The key simplification: $\frac{1}{2}(\frac{2}{s-1} - \frac{2}{s}) = \frac{1}{s-1} - \frac{1}{s} = \frac{1}{s(s-1)}$.
[/guided]
[/step]
[step:Verify the integral converges for all $s \in \mathbb{C}$ and deduce meromorphic continuation]
The function $\psi(t) = \sum_{n=1}^\infty e^{-\pi n^2 t}$ decays exponentially as $t \to \infty$: for $t \geq 1$, $\psi(t) \leq \sum_{n=1}^\infty e^{-\pi n^2} \cdot e^{-\pi(n^2-1)(t-1)}$, and in particular $\psi(t) \leq C e^{-\pi t}$ for some constant $C > 0$ and all $t \geq 1$.
For any $s \in \mathbb{C}$, the integrand $\left(t^{s/2-1} + t^{(1-s)/2-1}\right)\psi(t)$ on $[1, \infty)$ is bounded by
\begin{align*}
\left|t^{s/2-1} + t^{(1-s)/2-1}\right| \cdot |\psi(t)| &\leq \left(t^{\operatorname{Re}(s)/2-1} + t^{(1-\operatorname{Re}(s))/2-1}\right) \cdot Ce^{-\pi t}
\end{align*}
for $t \geq 1$. Since $e^{-\pi t}$ decays faster than any polynomial grows, the integral
\begin{align*}
\int_1^\infty \left(t^{s/2 - 1} + t^{(1-s)/2 - 1}\right)\psi(t)\, d\mathcal{L}^1(t)
\end{align*}
converges absolutely for every $s \in \mathbb{C}$, and defines an entire function of $s$ (by standard theorems on parameter-dependent integrals with dominated convergence — the dominated convergence theorem justifies differentiating under the integral sign).
Therefore
\begin{align*}
\xi(s) = \frac{1}{s(s-1)} + \int_1^\infty \left(t^{s/2 - 1} + t^{(1-s)/2 - 1}\right)\psi(t)\, d\mathcal{L}^1(t)
\end{align*}
provides the meromorphic continuation of $\xi(s) = \pi^{-s/2}\Gamma(s/2)\zeta(s)$ to all of $\mathbb{C}$, with poles only at $s = 0$ and $s = 1$ (from the term $1/(s(s-1))$).
Since $\Gamma(s/2)$ has a simple pole at $s = 0$ (with residue $2$), the factor $\pi^{-s/2}\Gamma(s/2)$ has a simple pole at $s = 0$. The pole of $\xi(s)$ at $s = 0$ is therefore accounted for by $\Gamma(s/2)$, and $\zeta(s)$ is regular at $s = 0$.
At $s = 1$: $\Gamma(1/2) = \sqrt{\pi}$ is finite, and $\pi^{-1/2}\Gamma(1/2) = \pi^{-1/2} \cdot \sqrt{\pi} = 1$. So $\xi(s)$ having a simple pole at $s = 1$ means $\zeta(s)$ has a simple pole at $s = 1$. The residue is
\begin{align*}
\operatorname{Res}_{s=1}\zeta(s) &= \lim_{s \to 1}(s-1)\zeta(s) = \lim_{s \to 1} \frac{(s-1)\xi(s)}{\pi^{-s/2}\Gamma(s/2)} = \frac{1}{\pi^{-1/2}\Gamma(1/2)} \cdot \lim_{s \to 1}\frac{(s-1)}{s(s-1)} \cdot s = \frac{1}{1} \cdot 1 = 1.
\end{align*}
More precisely: near $s = 1$, $\xi(s) = \frac{1}{s(s-1)} + (\text{holomorphic}) = \frac{-1}{s-1} \cdot \frac{1}{-s} + \cdots$. Actually $1/(s(s-1))$ has residue $-1$ at $s = 0$ and residue $1$ at $s = 1$. So $(s-1)\xi(s) \to 1$ as $s \to 1$. Since $\pi^{-s/2}\Gamma(s/2)\big|_{s=1} = \pi^{-1/2}\Gamma(1/2) = 1$, we get $\operatorname{Res}_{s=1}\zeta(s) = 1$.
[/step]
[step:Verify the functional equation from the symmetry of the representation]
The formula
\begin{align*}
\xi(s) &= \frac{1}{s(s-1)} + \int_1^\infty \left(t^{s/2 - 1} + t^{(1-s)/2 - 1}\right)\psi(t)\, d\mathcal{L}^1(t)
\end{align*}
is manifestly invariant under the substitution $s \mapsto 1 - s$: replacing $s$ by $1-s$ sends $s/2 - 1$ to $(1-s)/2 - 1$ and $(1-s)/2 - 1$ to $s/2 - 1$, so the integrand is unchanged. The rational term $1/(s(s-1)) = 1/((1-s)((1-s)-1))$ is also unchanged under $s \mapsto 1-s$ (since $s(s-1) = (1-s)(1-(1-s)) = (1-s)(-s) = s(s-1)$... more carefully: $(1-s)((1-s)-1) = (1-s)(-s) = s(s-1)$). Yes, $s(s-1) = (1-s)(-(1-s)+1-1) = (1-s)(-s)$... Let us verify: replacing $s$ by $1-s$ in $s(s-1)$ gives $(1-s)(1-s-1) = (1-s)(-s) = s(s-1)$. Correct.
Therefore $\xi(s) = \xi(1-s)$, which is precisely the functional equation
\begin{align*}
\pi^{-s/2}\,\Gamma\!\left(\frac{s}{2}\right)\zeta(s) &= \pi^{-(1-s)/2}\,\Gamma\!\left(\frac{1-s}{2}\right)\zeta(1-s).
\end{align*}
[guided]
The beauty of this proof is that the functional equation becomes a manifest symmetry of the final integral representation. The theta functional equation $\theta(1/t) = \sqrt{t}\,\theta(t)$ is the source of this symmetry: it is what allowed us to transform the integral over $(0,1)$ into an integral over $(1,\infty)$ with the exponent $s/2$ replaced by $(1-s)/2$.
Why is $s(s-1)$ invariant under $s \mapsto 1-s$? Because $s(s-1) = -(s)(1-s)$ and $(1-s)((1-s)-1) = (1-s)(-s) = -(1-s)(s)$, which equals $-(s)(1-s) = s(s-1)$. Both expressions equal $s^2 - s$, and substituting $s \to 1-s$ gives $(1-s)^2 - (1-s) = 1 - 2s + s^2 - 1 + s = s^2 - s$. So indeed the rational part is invariant.
The integrand $t^{s/2-1} + t^{(1-s)/2-1}$ is a sum of two terms that are exchanged by $s \mapsto 1-s$, so their sum is invariant. The function $\psi(t)$ and the domain $[1,\infty)$ are independent of $s$. Thus the entire expression is symmetric under $s \mapsto 1-s$, giving $\xi(s) = \xi(1-s)$.
[/guided]
[/step]
