[proofplan]
We first show that a weight two cusp form produces a $\Gamma_0(N)$-invariant holomorphic differential on the upper half-plane, hence a holomorphic differential on the non-cuspidal quotient. We then check the local behaviour at every cusp using the cusp parameter $q = e^{2\pi i z/h}$, where $h$ is the cusp width; the vanishing of the constant term is exactly what makes the differential holomorphic at $q=0$. Conversely, every holomorphic differential on $X_0(N)$ pulls back to a $\Gamma_0(N)$-invariant differential on $\mathcal{H}$, which has the form $2\pi i f(z)\,dz$ for a holomorphic weight two modular form, and holomorphy at the cusps forces $f$ to be cuspidal.
[/proofplan]
[step:Construct an invariant differential from a cusp form]
Let $f \in S_2(\Gamma_0(N))$. Define the holomorphic differential
\begin{align*}
\widetilde{\omega}_f := 2\pi i\, f(z)\,dz
\end{align*}
on $\mathcal{H}$. For
\begin{align*}
\gamma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in \Gamma_0(N),
\end{align*}
the action is given by
\begin{align*}
\gamma z = \frac{az+b}{cz+d},
\end{align*}
and
\begin{align*}
d(\gamma z) = \frac{1}{(cz+d)^2}\,dz.
\end{align*}
Since $f$ has weight $2$ for $\Gamma_0(N)$,
\begin{align*}
f(\gamma z) = (cz+d)^2 f(z).
\end{align*}
Therefore
\begin{align*}
\gamma^*\widetilde{\omega}_f
&= 2\pi i\, f(\gamma z)\,d(\gamma z) \\
&= 2\pi i\, (cz+d)^2 f(z)\frac{1}{(cz+d)^2}\,dz \\
&= 2\pi i\, f(z)\,dz \\
&= \widetilde{\omega}_f.
\end{align*}
Thus $\widetilde{\omega}_f$ is $\Gamma_0(N)$-invariant.
[guided]
The reason weight $2$ is the correct weight is visible in the differential transformation formula. If
\begin{align*}
\gamma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in \Gamma_0(N),
\end{align*}
then $\gamma$ acts on $\mathcal{H}$ by
\begin{align*}
z \mapsto \gamma z = \frac{az+b}{cz+d}.
\end{align*}
Differentiating gives
\begin{align*}
d(\gamma z) = \frac{1}{(cz+d)^2}\,dz.
\end{align*}
A modular form $f \in S_2(\Gamma_0(N))$ satisfies
\begin{align*}
f(\gamma z) = (cz+d)^2 f(z).
\end{align*}
The two factors cancel:
\begin{align*}
\gamma^*(2\pi i\, f(z)\,dz)
&= 2\pi i\, f(\gamma z)\,d(\gamma z) \\
&= 2\pi i\, (cz+d)^2 f(z)\frac{1}{(cz+d)^2}\,dz \\
&= 2\pi i\, f(z)\,dz.
\end{align*}
Thus the differential
\begin{align*}
\widetilde{\omega}_f := 2\pi i\, f(z)\,dz
\end{align*}
is invariant under every element of $\Gamma_0(N)$. This invariance is exactly the condition needed for the differential on $\mathcal{H}$ to descend to the quotient away from the cusps.
[/guided]
[/step]
[step:Descend the invariant differential across the non-cuspidal quotient]
Let $Y_0(N) := \Gamma_0(N)\backslash \mathcal{H}$, and let
\begin{align*}
\pi: \mathcal{H} \to Y_0(N)
\end{align*}
be the quotient map. Since $\widetilde{\omega}_f$ is $\Gamma_0(N)$-invariant, it defines a holomorphic differential on the smooth locus of $Y_0(N)$. At a point with finite stabilizer, we use the standard local normal form for a finite cyclic group of holomorphic automorphisms of a Riemann surface: there is a holomorphic coordinate $w$ on $\mathcal{H}$ centered at a lift of the point in which the stabilizer has order $e$ and acts by $w \mapsto \zeta w$, where $\zeta \in \mathbb{C}$ satisfies $\zeta^e = 1$ and has order $e$. Write the invariant local differential as
\begin{align*}
h(w)\,dw,
\end{align*}
where $h$ is holomorphic near $0$. If
\begin{align*}
h(w)=\sum_{m=0}^{\infty} b_m w^m,
\end{align*}
then invariance under $w \mapsto \zeta w$ gives
\begin{align*}
h(\zeta w)\zeta\,dw = h(w)\,dw,
\end{align*}
so $b_m(\zeta^{m+1}-1)=0$ for every $m \ge 0$. Hence only powers $m=re-1$ can occur. With the quotient coordinate $t=w^e$, we obtain
\begin{align*}
h(w)\,dw
= \sum_{r\ge 1} b_{re-1} w^{re-1}\,dw
= \frac{1}{e}\sum_{r\ge 1} b_{re-1} t^{r-1}\,dt,
\end{align*}
which is holomorphic in $t$. Therefore $\widetilde{\omega}_f$ descends to a holomorphic differential on $Y_0(N)$.
[/step]
[step:Verify holomorphy at every cusp]
Let $\mathfrak{a}$ be a cusp of $\Gamma_0(N)$. Choose
\begin{align*}
\sigma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z})
\end{align*}
with $\sigma\infty = \mathfrak{a}$, and let $h_{\mathfrak{a}}\in \mathbb{N}$ be the width of the cusp, so that the local parameter at $\mathfrak{a}$ is
\begin{align*}
q_{\mathfrak{a}} := e^{2\pi i z/h_{\mathfrak{a}}}.
\end{align*}
Define the weight-two slash transform
\begin{align*}
f|_2\sigma: \mathcal{H} &\to \mathbb{C} \\
z &\mapsto (cz+d)^{-2} f(\sigma z).
\end{align*}
The transformed cusp form $f|_2\sigma$ has Fourier expansion
\begin{align*}
(f|_2\sigma)(z) = \sum_{n=1}^{\infty} a_n q_{\mathfrak{a}}^n
\end{align*}
near $q_{\mathfrak{a}}=0$, because $f$ vanishes at the cusp $\mathfrak{a}$. Since
\begin{align*}
dq_{\mathfrak{a}} = \frac{2\pi i}{h_{\mathfrak{a}}} q_{\mathfrak{a}}\,dz,
\end{align*}
we have
\begin{align*}
2\pi i\,(f|_2\sigma)(z)\,dz
&= h_{\mathfrak{a}} (f|_2\sigma)(z)\frac{dq_{\mathfrak{a}}}{q_{\mathfrak{a}}} \\
&= h_{\mathfrak{a}}\sum_{n=1}^{\infty} a_n q_{\mathfrak{a}}^{n-1}\,dq_{\mathfrak{a}}.
\end{align*}
The final expression is holomorphic at $q_{\mathfrak{a}}=0$. Thus the descended differential extends holomorphically across every cusp, giving an element
\begin{align*}
\omega_f \in H^0(X_0(N), \Omega^1_{X_0(N)}).
\end{align*}
[guided]
We now check the only possible singularities of the descended differential: the cusps. Fix a cusp $\mathfrak{a}$. Choose
\begin{align*}
\sigma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z})
\end{align*}
carrying $\infty$ to $\mathfrak{a}$, and let $h_{\mathfrak{a}}$ be the width of the cusp. In this coordinate, the quotient near the cusp is described by the parameter
\begin{align*}
q_{\mathfrak{a}} := e^{2\pi i z/h_{\mathfrak{a}}}.
\end{align*}
Define the weight-two slash transform by
\begin{align*}
f|_2\sigma: \mathcal{H} &\to \mathbb{C} \\
z &\mapsto (cz+d)^{-2} f(\sigma z).
\end{align*}
The cusp condition says precisely that the Fourier expansion of $f$ at $\mathfrak{a}$ has no constant term. Equivalently,
\begin{align*}
(f|_2\sigma)(z) = \sum_{n=1}^{\infty} a_n q_{\mathfrak{a}}^n.
\end{align*}
The differential relation between $z$ and the cusp parameter is
\begin{align*}
dq_{\mathfrak{a}} = \frac{2\pi i}{h_{\mathfrak{a}}}q_{\mathfrak{a}}\,dz.
\end{align*}
Solving for $2\pi i\,dz$ gives
\begin{align*}
2\pi i\,dz = h_{\mathfrak{a}}\frac{dq_{\mathfrak{a}}}{q_{\mathfrak{a}}}.
\end{align*}
Therefore
\begin{align*}
2\pi i\,(f|_2\sigma)(z)\,dz
&= h_{\mathfrak{a}}(f|_2\sigma)(z)\frac{dq_{\mathfrak{a}}}{q_{\mathfrak{a}}} \\
&= h_{\mathfrak{a}}\left(\sum_{n=1}^{\infty} a_n q_{\mathfrak{a}}^n\right)\frac{dq_{\mathfrak{a}}}{q_{\mathfrak{a}}} \\
&= h_{\mathfrak{a}}\sum_{n=1}^{\infty} a_n q_{\mathfrak{a}}^{n-1}\,dq_{\mathfrak{a}}.
\end{align*}
This is a convergent [power series](/page/Power%20Series) in $q_{\mathfrak{a}}$ multiplied by $dq_{\mathfrak{a}}$, so it is holomorphic at $q_{\mathfrak{a}}=0$. The factor $2\pi i$ is exactly what converts $dz$ into $dq/q$; this is why the normalization matches Fourier expansions in the variable $q$.
[/guided]
[/step]
[step:Recover a weight two cusp form from a holomorphic differential]
Let
\begin{align*}
\omega \in H^0(X_0(N), \Omega^1_{X_0(N)}).
\end{align*}
Pull it back to $\mathcal{H}$ along $\pi$:
\begin{align*}
\widetilde{\omega} := \pi^*\omega.
\end{align*}
Since $dz$ is a nowhere-vanishing holomorphic differential on $\mathcal{H}$, there is a unique [holomorphic function](/page/Holomorphic%20Function)
\begin{align*}
f: \mathcal{H} \to \mathbb{C}
\end{align*}
such that
\begin{align*}
\widetilde{\omega} = 2\pi i\, f(z)\,dz.
\end{align*}
For every $\gamma \in \Gamma_0(N)$, the equality $\pi\circ\gamma=\pi$ gives
\begin{align*}
\gamma^*\widetilde{\omega} = \widetilde{\omega}.
\end{align*}
Writing $\gamma=\begin{pmatrix}a&b\\ c&d\end{pmatrix}$, this means
\begin{align*}
2\pi i\, f(\gamma z)\frac{1}{(cz+d)^2}\,dz
= 2\pi i\, f(z)\,dz.
\end{align*}
Hence
\begin{align*}
f(\gamma z) = (cz+d)^2 f(z),
\end{align*}
so $f$ has weight $2$ for $\Gamma_0(N)$.
It remains to check cuspidality. At a cusp $\mathfrak{a}$, choose
\begin{align*}
\sigma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z})
\end{align*}
with $\sigma\infty=\mathfrak{a}$. Let $h_{\mathfrak{a}}\in\mathbb{N}$ be the cusp width, and define the local parameter $q_{\mathfrak{a}}:=e^{2\pi i z/h_{\mathfrak{a}}}$. Define
\begin{align*}
f|_2\sigma: \mathcal{H} &\to \mathbb{C} \\
z &\mapsto (cz+d)^{-2} f(\sigma z).
\end{align*}
Holomorphy of $\omega$ gives a convergent expansion
\begin{align*}
\omega = \left(\sum_{m=0}^{\infty} b_m q_{\mathfrak{a}}^m\right)dq_{\mathfrak{a}}.
\end{align*}
Pulling back to the upper half-plane gives
\begin{align*}
2\pi i\,(f|_2\sigma)(z)\,dz
= \left(\sum_{m=0}^{\infty} b_m q_{\mathfrak{a}}^m\right)dq_{\mathfrak{a}}
= \left(\sum_{m=0}^{\infty} b_m q_{\mathfrak{a}}^m\right)\frac{2\pi i}{h_{\mathfrak{a}}}q_{\mathfrak{a}}\,dz.
\end{align*}
Thus
\begin{align*}
(f|_2\sigma)(z)
= \frac{1}{h_{\mathfrak{a}}}\sum_{m=0}^{\infty} b_m q_{\mathfrak{a}}^{m+1}.
\end{align*}
The constant term is zero at every cusp, so $f \in S_2(\Gamma_0(N))$.
[guided]
Start with a holomorphic differential
\begin{align*}
\omega \in H^0(X_0(N), \Omega^1_{X_0(N)}).
\end{align*}
Pull it back to the upper half-plane by the quotient map $\pi: \mathcal{H}\to Y_0(N)$ and write
\begin{align*}
\widetilde{\omega} := \pi^*\omega.
\end{align*}
Because $dz$ is a nowhere-vanishing holomorphic differential on $\mathcal{H}$, there is a unique holomorphic function
\begin{align*}
f: \mathcal{H} &\to \mathbb{C}
\end{align*}
such that
\begin{align*}
\widetilde{\omega}=2\pi i\,f(z)\,dz.
\end{align*}
The pullback is invariant under $\Gamma_0(N)$, since $\pi\circ\gamma=\pi$ for every $\gamma\in\Gamma_0(N)$. Thus
\begin{align*}
\gamma^*\widetilde{\omega}=\widetilde{\omega}.
\end{align*}
Writing
\begin{align*}
\gamma=\begin{pmatrix}a&b\\ c&d\end{pmatrix},
\end{align*}
we compute
\begin{align*}
2\pi i\, f(\gamma z)\frac{1}{(cz+d)^2}\,dz
=2\pi i\,f(z)\,dz.
\end{align*}
Cancelling the nowhere-zero factor $2\pi i\,dz$ gives
\begin{align*}
f(\gamma z)=(cz+d)^2f(z),
\end{align*}
so $f$ is a weight two modular form for $\Gamma_0(N)$.
It remains to see that $f$ vanishes at every cusp. Fix a cusp $\mathfrak{a}$, choose
\begin{align*}
\sigma=\begin{pmatrix}a&b\\ c&d\end{pmatrix}\in SL_2(\mathbb{Z})
\end{align*}
with $\sigma\infty=\mathfrak{a}$, and let $h_{\mathfrak{a}}\in\mathbb{N}$ be the cusp width. Define
\begin{align*}
q_{\mathfrak{a}}:=e^{2\pi i z/h_{\mathfrak{a}}}
\end{align*}
and define the weight-two slash transform by
\begin{align*}
f|_2\sigma: \mathcal{H} &\to \mathbb{C} \\
z &\mapsto (cz+d)^{-2}f(\sigma z).
\end{align*}
Since $\omega$ is holomorphic at the cusp, it has a convergent local expansion
\begin{align*}
\omega=\left(\sum_{m=0}^{\infty}b_m q_{\mathfrak{a}}^m\right)dq_{\mathfrak{a}}.
\end{align*}
The differential relation
\begin{align*}
dq_{\mathfrak{a}}=\frac{2\pi i}{h_{\mathfrak{a}}}q_{\mathfrak{a}}\,dz
\end{align*}
therefore gives
\begin{align*}
2\pi i\,(f|_2\sigma)(z)\,dz
&=\left(\sum_{m=0}^{\infty}b_m q_{\mathfrak{a}}^m\right)dq_{\mathfrak{a}} \\
&=\left(\sum_{m=0}^{\infty}b_m q_{\mathfrak{a}}^m\right)\frac{2\pi i}{h_{\mathfrak{a}}}q_{\mathfrak{a}}\,dz.
\end{align*}
Cancelling $2\pi i\,dz$ yields
\begin{align*}
(f|_2\sigma)(z)=\frac{1}{h_{\mathfrak{a}}}\sum_{m=0}^{\infty}b_m q_{\mathfrak{a}}^{m+1}.
\end{align*}
This expansion has no constant term. Since the cusp $\mathfrak{a}$ was arbitrary, $f$ vanishes at every cusp, and hence $f\in S_2(\Gamma_0(N))$.
[/guided]
[/step]
[step:Check that the two constructions are inverse and complex-linear]
The construction $f \mapsto \omega_f$ is complex-linear because
\begin{align*}
\omega_{\alpha f+\beta g}
= 2\pi i\,(\alpha f+\beta g)(z)\,dz
= \alpha\,\omega_f+\beta\,\omega_g
\end{align*}
for all $\alpha,\beta\in\mathbb{C}$ and all $f,g\in S_2(\Gamma_0(N))$. Starting with $f$, descending $2\pi i f(z)\,dz$ and then pulling back to $\mathcal{H}$ recovers $2\pi i f(z)\,dz$, hence recovers $f$. Starting with $\omega$, pulling back to $2\pi i f(z)\,dz$ and then descending gives the unique differential on $X_0(N)$ whose pullback is $\pi^*\omega$, namely $\omega$. Therefore
\begin{align*}
S_2(\Gamma_0(N)) \longrightarrow H^0(X_0(N), \Omega^1_{X_0(N)}), \qquad f \longmapsto \omega_f,
\end{align*}
is a natural complex-linear isomorphism.
[guided]
The assignment is complex-linear because for all $\alpha,\beta\in\mathbb{C}$ and all $f,g\in S_2(\Gamma_0(N))$,
\begin{align*}
\omega_{\alpha f+\beta g}
&=2\pi i\,(\alpha f+\beta g)(z)\,dz \\
&=\alpha\,2\pi i\,f(z)\,dz+\beta\,2\pi i\,g(z)\,dz \\
&=\alpha\,\omega_f+\beta\,\omega_g.
\end{align*}
Now check that the two constructions undo each other. If we start with a cusp form $f$, construct $2\pi i f(z)\,dz$ on $\mathcal{H}$, descend it to $X_0(N)$, and pull it back again along $\pi$, the pullback is the original invariant differential $2\pi i f(z)\,dz$. Since $dz$ is nowhere-vanishing on $\mathcal{H}$, this recovers the original function $f$.
Conversely, if we start with a holomorphic differential $\omega$ on $X_0(N)$, pull it back to $\widetilde{\omega}=2\pi i f(z)\,dz$, and then descend that invariant differential, the descended differential has pullback $\pi^*\omega$. The quotient construction gives uniqueness of the descended differential with this pullback, so the result is exactly $\omega$. Hence the map
\begin{align*}
S_2(\Gamma_0(N)) \longrightarrow H^0(X_0(N), \Omega^1_{X_0(N)}), \qquad f \longmapsto \omega_f,
\end{align*}
is a complex-linear bijection, and this is the claimed natural complex-linear isomorphism.
[/guided]
[/step]