[proofplan]
The proof realizes the completed $L$-function as the Mellin transform of $f$ along the positive imaginary axis. The modular relation under $z \mapsto -1/z$ converts the part of the Mellin integral over $(0,1)$ into another integral over $(1,\infty)$, producing an entire expression symmetric under $s \mapsto k-s$. Finally, the Fourier expansion justifies the equality between this Mellin transform and $(2\pi)^{-s}\Gamma(s)L(f,s)$ in a right half-plane, so the entire expression is the desired [analytic continuation](/page/Analytic%20Continuation).
[/proofplan]
[step:Reduce to even weight or the zero form]
If $k$ is odd, then $f=0$. Indeed, the matrix $-I \in SL_2(\mathbb{Z})$ fixes every point of $\mathbb{H}$, while the weight-$k$ modularity relation gives
\begin{align*}
f(z) = f((-I)z) = (-1)^k f(z)
\end{align*}
for every $z \in \mathbb{H}$. If $k$ is odd, this implies $f(z)=0$ for all $z \in \mathbb{H}$, and the theorem is immediate.
We therefore assume for the rest of the proof that either $f=0$ or $k$ is even. In the nonzero case, $i^{-k}=i^k$ and $i^{2k}=1$.
[/step]
[step:Control the Fourier coefficients and justify the Dirichlet series]
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. Let $a_n \in \mathbb{C}$ denote the $n$th Fourier coefficient of $f$ at infinity. We first record a polynomial bound for $(a_n)_{n \geq 1}$. Define
\begin{align*}
\Phi: \mathbb{H} &\to [0,\infty) \\
z &\mapsto (\operatorname{Im} z)^{k/2}|f(z)|.
\end{align*}
For $\gamma=\begin{pmatrix}a&b\\ c&d\end{pmatrix} \in SL_2(\mathbb{Z})$, the identities
\begin{align*}
\operatorname{Im}(\gamma z)=\frac{\operatorname{Im} z}{|cz+d|^2},
\qquad
f(\gamma z)=(cz+d)^k f(z)
\end{align*}
show that $\Phi(\gamma z)=\Phi(z)$. Since $f$ is cuspidal at infinity, $\Phi(z)\to 0$ as $\operatorname{Im} z \to \infty$ in the standard fundamental domain. Choose $Y_0>0$ such that $\Phi(z)\leq 1$ on the part of the standard fundamental domain with $\operatorname{Im} z\geq Y_0$. The truncated part with $\operatorname{Im} z\leq Y_0$ is compact after taking the usual closed fundamental-domain representative, and $\Phi$ is continuous there, so $\Phi$ has a finite maximum $C_0\geq 0$ on that truncated part. Since every $z\in\mathbb{H}$ is $SL_2(\mathbb{Z})$-equivalent to a point in the standard fundamental domain and $\Phi$ is invariant under this action, defining $C_f:=\max\{1,C_0\}$ gives
\begin{align*}
|f(x+iy)| \leq C_f y^{-k/2}
\end{align*}
for all $x \in \mathbb{R}$ and $y>0$.
For $n \geq 1$, the Fourier coefficient may be recovered at height $y>0$ by
\begin{align*}
a_n = e^{2\pi n y}\int_0^1 f(x+iy)e^{-2\pi i n x}\,d\mathcal{L}^1(x).
\end{align*}
Taking $y=1/n$ gives
\begin{align*}
|a_n|
&\leq e^{2\pi} \int_0^1 |f(x+i/n)|\,d\mathcal{L}^1(x) \\
&\leq e^{2\pi} C_f n^{k/2}.
\end{align*}
Therefore the series
\begin{align*}
L(f,s)=\sum_{n=1}^{\infty}a_n n^{-s}
\end{align*}
converges absolutely whenever $\operatorname{Re}(s)>1+k/2$.
[/step]
[step:Identify the completed $L$-function with a Mellin transform]
Define the Mellin transform on the right half-plane by
\begin{align*}
M_f:\{s \in \mathbb{C}: \operatorname{Re}(s)>1+k/2\} &\to \mathbb{C} \\
s &\mapsto \int_0^\infty f(iy)y^{s-1}\,d\mathcal{L}^1(y).
\end{align*}
Using the Fourier expansion
\begin{align*}
f(iy)=\sum_{n=1}^{\infty}a_n e^{-2\pi n y},
\end{align*}
we justify termwise integration by Tonelli's theorem applied to the non-negative summands $|a_n|e^{-2\pi n y}y^{\operatorname{Re}(s)-1}$. The coefficient bound from the previous step gives
\begin{align*}
\sum_{n=1}^{\infty}|a_n|\int_0^\infty e^{-2\pi n y}y^{\operatorname{Re}(s)-1}\,d\mathcal{L}^1(y)
&\leq e^{2\pi}C_f\sum_{n=1}^{\infty} n^{k/2}\int_0^\infty e^{-2\pi n y}y^{\operatorname{Re}(s)-1}\,d\mathcal{L}^1(y) \\
&= e^{2\pi}C_f(2\pi)^{-\operatorname{Re}(s)}\Gamma(\operatorname{Re}(s))
\sum_{n=1}^{\infty} n^{k/2-\operatorname{Re}(s)}.
\end{align*}
The final series converges because $\operatorname{Re}(s)>1+k/2$. Hence [Fubini's theorem](/theorems/2961) applies to the absolutely integrable series and yields
\begin{align*}
M_f(s)
&=\sum_{n=1}^{\infty}a_n\int_0^\infty e^{-2\pi n y}y^{s-1}\,d\mathcal{L}^1(y).
\end{align*}
In the inner integral, use the substitution $u=2\pi n y$, so that $d\mathcal{L}^1(u)=2\pi n\,d\mathcal{L}^1(y)$ and $y=u/(2\pi n)$. Then
\begin{align*}
\int_0^\infty e^{-2\pi n y}y^{s-1}\,d\mathcal{L}^1(y)
&=(2\pi n)^{-s}\int_0^\infty e^{-u}u^{s-1}\,d\mathcal{L}^1(u) \\
&=(2\pi n)^{-s}\Gamma(s).
\end{align*}
Therefore
\begin{align*}
M_f(s)
&=(2\pi)^{-s}\Gamma(s)\sum_{n=1}^{\infty}a_n n^{-s} \\
&=\Lambda(f,s)
\end{align*}
for $\operatorname{Re}(s)>1+k/2$.
[/step]
[step:Rewrite the Mellin transform as an entire symmetric integral]
Split the Mellin integral at $1$:
\begin{align*}
M_f(s)
=
\int_1^\infty f(iy)y^{s-1}\,d\mathcal{L}^1(y)
+
\int_0^1 f(iy)y^{s-1}\,d\mathcal{L}^1(y).
\end{align*}
In the second integral, use the substitution $t=1/y$, so that $y=1/t$, $d\mathcal{L}^1(y)=t^{-2}\,d\mathcal{L}^1(t)$ with the orientation reversed, and the interval $(0,1)$ becomes $(1,\infty)$. The modular relation for $S=\begin{pmatrix}0&-1\\1&0\end{pmatrix}$ gives
\begin{align*}
f(i/t)=f(S(it))=(it)^k f(it).
\end{align*}
Equivalently,
\begin{align*}
f(iy)=i^{-k}y^{-k}f(i/y).
\end{align*}
Since $k$ is even in the nonzero case, $i^{-k}=i^k$. Hence
\begin{align*}
\int_0^1 f(iy)y^{s-1}\,d\mathcal{L}^1(y)
&=i^k\int_1^\infty f(it)t^{k-s-1}\,d\mathcal{L}^1(t).
\end{align*}
Renaming $t$ as $y$, define
\begin{align*}
F_f:\mathbb{C} &\to \mathbb{C} \\
s &\mapsto \int_1^\infty f(iy)\left(y^{s-1}+i^k y^{k-s-1}\right)\,d\mathcal{L}^1(y).
\end{align*}
We now justify the analytic dependence on $s$. From the Fourier expansion and the coefficient bound already proved, for $y\geq 1$ we have
\begin{align*}
|f(iy)|
&\leq \sum_{n=1}^{\infty}|a_n|e^{-2\pi n y} \\
&\leq e^{2\pi}C_f\sum_{n=1}^{\infty} n^{k/2}e^{-2\pi n y} \\
&\leq e^{2\pi}C_f e^{-\pi y}\sum_{n=1}^{\infty} n^{k/2}e^{-\pi n}.
\end{align*}
Define the finite constant
\begin{align*}
C_1:=e^{2\pi}C_f\sum_{n=1}^{\infty} n^{k/2}e^{-\pi n}.
\end{align*}
Then $|f(iy)|\leq C_1e^{-\pi y}$ for all $y\geq 1$. Let $K\subset\mathbb{C}$ be compact, and choose $A_K>0$ such that $|\operatorname{Re}(s)|\leq A_K$ for all $s\in K$. For $y\geq 1$ and $s\in K$,
\begin{align*}
\left|f(iy)\left(y^{s-1}+i^k y^{k-s-1}\right)\right|
&\leq C_1e^{-\pi y}\left(y^{A_K-1}+y^{k+A_K-1}\right).
\end{align*}
The dominating function on the right is integrable on $(1,\infty)$ with respect to $\mathcal{L}^1$. Since for each fixed $y\geq 1$ the map $s\mapsto f(iy)(y^{s-1}+i^k y^{k-s-1})$ is entire, the holomorphic parameter integral theorem implies that $F_f$ is entire. Since $F_f(s)=M_f(s)=\Lambda(f,s)$ on the half-plane $\operatorname{Re}(s)>1+k/2$, $F_f$ is the analytic continuation of $\Lambda(f,s)$ to an entire function.
[/step]
[step:Derive the functional equation from the symmetric integral]
For every $s \in \mathbb{C}$, the entire continuation is
\begin{align*}
\Lambda(f,s)
=
F_f(s)
=
\int_1^\infty f(iy)\left(y^{s-1}+i^k y^{k-s-1}\right)\,d\mathcal{L}^1(y).
\end{align*}
Similarly,
\begin{align*}
F_f(k-s)
&=
\int_1^\infty f(iy)\left(y^{k-s-1}+i^k y^{s-1}\right)\,d\mathcal{L}^1(y).
\end{align*}
Multiplying by $i^k$ and using $i^{2k}=1$ in the nonzero case gives
\begin{align*}
i^k F_f(k-s)
&=
\int_1^\infty f(iy)\left(i^k y^{k-s-1}+i^{2k}y^{s-1}\right)\,d\mathcal{L}^1(y) \\
&=
\int_1^\infty f(iy)\left(y^{s-1}+i^k y^{k-s-1}\right)\,d\mathcal{L}^1(y) \\
&=F_f(s).
\end{align*}
Thus
\begin{align*}
\Lambda(f,s)=i^k\Lambda(f,k-s)
\end{align*}
for all $s \in \mathbb{C}$, completing the proof.
[/step]