[proofplan]
We restrict the Fourier expansion of the cusp form to the positive imaginary axis, where $q=e^{-2\pi y}$. The cusp condition gives exponential decay as $y\to\infty$, while modularity under $z\mapsto -1/z$ converts the behavior near $y=0$ into exponential decay at infinity. For $\operatorname{Re}(s)$ sufficiently large, absolute convergence of the Dirichlet series and the Mellin integral allow termwise integration, and each exponential term evaluates to the Gamma integral after the substitution $t=2\pi n y$.
[/proofplan]
[step:Restrict the Fourier expansion to the imaginary axis]
For $y \in (0,\infty)$, define the restriction
\begin{align*}
F:(0,\infty)&\to \mathbb{C}\\
y&\mapsto f(iy).
\end{align*}
Since $e^{2\pi i n(iy)}=e^{-2\pi n y}$, the Fourier expansion of $f$ gives
\begin{align*}
F(y)=f(iy)=\sum_{n=1}^{\infty} a_n e^{-2\pi n y}.
\end{align*}
The convergence is locally uniform for $y>0$, because the Fourier expansion of a holomorphic modular form converges normally on sets bounded away from the cusp.
[guided]
We first translate the modular-form Fourier expansion into a real-variable formula along the positive imaginary axis. Define
\begin{align*}
F:(0,\infty)&\to \mathbb{C}\\
y&\mapsto f(iy).
\end{align*}
The variable $y$ is positive, so $iy \in \mathbb{H}$. In the Fourier expansion
\begin{align*}
f(z)=\sum_{n=1}^{\infty}a_n e^{2\pi i n z},
\end{align*}
we substitute $z=iy$. The exponential factor becomes
\begin{align*}
e^{2\pi i n(iy)}=e^{-2\pi n y}.
\end{align*}
Therefore
\begin{align*}
F(y)=f(iy)=\sum_{n=1}^{\infty}a_n e^{-2\pi n y}.
\end{align*}
The convergence is locally uniform on $(0,\infty)$ because every compact subinterval of $(0,\infty)$ stays a positive distance away from the cusp $y=\infty$ in the $q$-coordinate and the Fourier expansion converges normally there.
[/guided]
[/step]
[step:Record the convergence needed for the Mellin integral]
Because $f$ is cuspidal at $\infty$, there exist constants $C_\infty>0$ and $\eta_\infty>0$ such that
\begin{align*}
|f(iy)|\le C_\infty e^{-\eta_\infty y}
\end{align*}
for all $y\ge 1$. Since $f$ has weight $k$ for $SL_2(\mathbb{Z})$, the modular transformation $S:z\mapsto -1/z$ gives
\begin{align*}
f(i/y)=(iy)^k f(iy)
\end{align*}
for $0<y\le 1$. Hence
\begin{align*}
|f(iy)|=y^{-k}|f(i/y)|.
\end{align*}
Applying the cusp decay estimate to $i/y$ gives constants $C_0>0$ and $\eta_0>0$ such that
\begin{align*}
|f(iy)|\le C_0 y^{-k} e^{-\eta_0/y}
\end{align*}
for $0<y\le 1$.
Thus, for every $s\in\mathbb{C}$, the function
\begin{align*}
M_s:(0,\infty)&\to \mathbb{C}\\
y&\mapsto f(iy)y^s y^{-1}
\end{align*}
is absolutely integrable with respect to $\mathcal{L}^1$ on both $(0,1]$ and $[1,\infty)$. In particular, the Mellin integral
\begin{align*}
\int_0^\infty f(iy)y^s\,\frac{d\mathcal{L}^1(y)}{y}
\end{align*}
is well-defined for every $s\in\mathbb{C}$.
[guided]
The Mellin integral has two possible problem points: $y=\infty$ and $y=0$. At infinity, the cusp condition says that the constant term in the Fourier expansion vanishes, so the first nonzero exponential term already decays exponentially. Therefore there are constants $C_\infty>0$ and $\eta_\infty>0$ such that
\begin{align*}
|f(iy)|\le C_\infty e^{-\eta_\infty y}
\end{align*}
for every $y\ge 1$.
Near $y=0$, we use modularity to move the point $iy$ back to height $1/y$, which is large. The matrix
\begin{align*}
S=\begin{pmatrix}0&-1\\1&0\end{pmatrix}
\end{align*}
acts by $S z=-1/z$. Since $f$ has weight $k$ for $SL_2(\mathbb{Z})$,
\begin{align*}
f(Sz)=z^k f(z).
\end{align*}
Substituting $z=iy$ gives $S(iy)=i/y$, and hence
\begin{align*}
f(i/y)=(iy)^k f(iy).
\end{align*}
Taking absolute values yields
\begin{align*}
|f(iy)|=y^{-k}|f(i/y)|.
\end{align*}
Now $1/y\ge 1$ when $0<y\le 1$, so the cusp decay estimate applies to $f(i/y)$. Thus there are constants $C_0>0$ and $\eta_0>0$ such that
\begin{align*}
|f(iy)|\le C_0 y^{-k}e^{-\eta_0/y}
\end{align*}
for every $0<y\le 1$.
Define
\begin{align*}
M_s:(0,\infty)&\to \mathbb{C}\\
y&\mapsto f(iy)y^s y^{-1}.
\end{align*}
On $[1,\infty)$, the absolute value of $M_s$ is bounded by an exponential times a power of $y$, which is integrable with respect to $\mathcal{L}^1$. On $(0,1]$, the factor $e^{-\eta_0/y}$ decays faster than every power of $y^{-1}$, so $y^{\operatorname{Re}(s)-k-1}e^{-\eta_0/y}$ is integrable with respect to $\mathcal{L}^1$. Hence
\begin{align*}
\int_0^\infty |f(iy)y^s|\,\frac{d\mathcal{L}^1(y)}{y}<\infty.
\end{align*}
This proves that the Mellin integral is well-defined.
[/guided]
[/step]
[step:Justify termwise integration in a right half-plane]
Choose a real number $\sigma_0$ such that the Dirichlet series
\begin{align*}
\sum_{n=1}^{\infty} |a_n|n^{-\sigma}
\end{align*}
converges for every $\sigma>\sigma_0$. Let $s\in\mathbb{C}$ satisfy $\operatorname{Re}(s)=\sigma>\max\{\sigma_0,0\}$. For each $n\in\mathbb{N}$, define
\begin{align*}
G_{n,s}:(0,\infty)&\to \mathbb{C}\\
y&\mapsto a_n e^{-2\pi n y}y^s y^{-1}.
\end{align*}
Then
\begin{align*}
\int_0^\infty |G_{n,s}(y)|\,d\mathcal{L}^1(y)
=
|a_n|\int_0^\infty e^{-2\pi n y}y^{\sigma}\,\frac{d\mathcal{L}^1(y)}{y}.
\end{align*}
By the Gamma integral computation in the next step, this equals
\begin{align*}
|a_n|(2\pi n)^{-\sigma}\Gamma(\sigma).
\end{align*}
Therefore
\begin{align*}
\sum_{n=1}^{\infty}\int_0^\infty |G_{n,s}(y)|\,d\mathcal{L}^1(y)
=
(2\pi)^{-\sigma}\Gamma(\sigma)\sum_{n=1}^{\infty}|a_n|n^{-\sigma}<\infty.
\end{align*}
By absolute summability in $L^1((0,\infty),\mathcal{L}^1)$, the series $\sum_n G_{n,s}$ may be integrated term by term. Hence
\begin{align*}
\int_0^\infty f(iy)y^s\,\frac{d\mathcal{L}^1(y)}{y}
=
\sum_{n=1}^{\infty}a_n\int_0^\infty e^{-2\pi n y}y^s\,\frac{d\mathcal{L}^1(y)}{y}.
\end{align*}
[guided]
We now explain why the integral of the [Fourier series](/page/Fourier%20Series) is the sum of the integrals. Choose $\sigma_0\in\mathbb{R}$ such that
\begin{align*}
\sum_{n=1}^{\infty}|a_n|n^{-\sigma}<\infty
\end{align*}
whenever $\sigma>\sigma_0$. Such a number exists because the Fourier coefficients of a holomorphic cusp form have at most polynomial growth. Fix $s\in\mathbb{C}$ and write $\sigma=\operatorname{Re}(s)$, assuming $\sigma>\max\{\sigma_0,0\}$.
For each $n\in\mathbb{N}$, define the measurable function
\begin{align*}
G_{n,s}:(0,\infty)&\to \mathbb{C}\\
y&\mapsto a_n e^{-2\pi n y}y^s y^{-1}.
\end{align*}
The absolute value is
\begin{align*}
|G_{n,s}(y)|=|a_n|e^{-2\pi n y}y^\sigma y^{-1}.
\end{align*}
Therefore
\begin{align*}
\int_0^\infty |G_{n,s}(y)|\,d\mathcal{L}^1(y)
=
|a_n|\int_0^\infty e^{-2\pi n y}y^\sigma\,\frac{d\mathcal{L}^1(y)}{y}.
\end{align*}
The integral on the right is the usual Gamma integral after the substitution $t=2\pi n y$, and equals $(2\pi n)^{-\sigma}\Gamma(\sigma)$. Thus
\begin{align*}
\sum_{n=1}^{\infty}\int_0^\infty |G_{n,s}(y)|\,d\mathcal{L}^1(y)
=
(2\pi)^{-\sigma}\Gamma(\sigma)\sum_{n=1}^{\infty}|a_n|n^{-\sigma}<\infty.
\end{align*}
This proves that $\sum_n G_{n,s}$ is absolutely summable in $L^1((0,\infty),\mathcal{L}^1)$. Consequently, the integral of the sum is the sum of the integrals. Since Step 1 identified
\begin{align*}
f(iy)y^s y^{-1}=\sum_{n=1}^{\infty}G_{n,s}(y),
\end{align*}
we obtain
\begin{align*}
\int_0^\infty f(iy)y^s\,\frac{d\mathcal{L}^1(y)}{y}
=
\sum_{n=1}^{\infty}a_n\int_0^\infty e^{-2\pi n y}y^s\,\frac{d\mathcal{L}^1(y)}{y}.
\end{align*}
[/guided]
[/step]
[step:Evaluate each exponential Mellin integral]
Fix $n\in\mathbb{N}$ and $s\in\mathbb{C}$ with $\operatorname{Re}(s)>0$. Define
\begin{align*}
I_n(s)=\int_0^\infty e^{-2\pi n y}y^s\,\frac{d\mathcal{L}^1(y)}{y}.
\end{align*}
Use the change of variables $t=2\pi n y$, so $y=t/(2\pi n)$ and $d\mathcal{L}^1(y)=(2\pi n)^{-1}d\mathcal{L}^1(t)$. The interval $(0,\infty)$ maps onto $(0,\infty)$. Hence
\begin{align*}
I_n(s)
&=\int_0^\infty e^{-t}\left(\frac{t}{2\pi n}\right)^s\frac{d\mathcal{L}^1(t)}{t}\\
&=(2\pi n)^{-s}\int_0^\infty e^{-t}t^s\,\frac{d\mathcal{L}^1(t)}{t}\\
&=(2\pi n)^{-s}\Gamma(s).
\end{align*}
[/step]
[step:Sum the evaluated terms to obtain the Mellin formula]
Substituting the computation of $I_n(s)$ into the termwise integral identity gives
\begin{align*}
\int_0^\infty f(iy)y^s\,\frac{d\mathcal{L}^1(y)}{y}
&=
\sum_{n=1}^{\infty}a_n(2\pi n)^{-s}\Gamma(s)\\
&=(2\pi)^{-s}\Gamma(s)\sum_{n=1}^{\infty}a_n n^{-s}\\
&=(2\pi)^{-s}\Gamma(s)L(f,s).
\end{align*}
This is the desired identity for every $s$ in the chosen right half-plane of absolute convergence, hence for $\operatorname{Re}(s)$ sufficiently large.
[/step]
