[step:Identify cusp forms with holomorphic differentials on the modular curve]
Let $\mathfrak{H}=\{z\in\mathbb{C}:\operatorname{Im}(z)>0\}$ denote the complex upper half-plane. Let
\begin{align*}
\Phi: S_2(\Gamma_0(N)) &\to H^0(X_0(N), \Omega^1) \\
f &\mapsto 2\pi i\, f(z)\, dz
\end{align*}
be the map in the statement. If $f \in S_2(\Gamma_0(N))$, then the weight $2$ transformation law gives
\begin{align*}
f(\delta z)\, d(\delta z) = f(z)\, dz
\end{align*}
for every $\delta \in \Gamma_0(N)$, because $d(\delta z) = (cz+d)^{-2}\,dz$ when $\delta=\begin{pmatrix}a&b\\ c&d\end{pmatrix}$. Hence $2\pi i f(z)\,dz$ descends from $\mathfrak{H}$ to a holomorphic differential on the open modular curve.
At a cusp $s$, choose a scaling matrix $\sigma_s\in SL_2(\mathbb{R})$ sending $\infty$ to $s$, and let $w_s\in\mathbb{N}$ denote the width of $s$ for $\Gamma_0(N)$. Define the local parameter
\begin{align*}
q_s: \sigma_s^{-1}\mathfrak{H}/\langle z\mapsto z+w_s\rangle &\to \mathbb{C} \\
z &\mapsto e^{2\pi i z/w_s}.
\end{align*}
Since $dq_s=(2\pi i/w_s)q_s\,dz$ in this cusp coordinate, and since cuspidality gives an expansion
\begin{align*}
f(\sigma_s z)(c_sz+d_s)^{-2}=\sum_{m=1}^{\infty} a_{s,m} q_s^m
\end{align*}
for the corresponding scaling matrix entries $c_s,d_s$, the descended differential has the form
\begin{align*}
2\pi i f(z)\,dz = w_s\sum_{m=1}^{\infty} a_{s,m} q_s^{m-1}\,dq_s
\end{align*}
in the parameter $q_s$. This is holomorphic at $s$ because the Fourier expansion has no constant term. Thus $\Phi(f)$ extends to $X_0(N)$. Conversely, every holomorphic differential on $X_0(N)$ restricts to a $\Gamma_0(N)$-invariant differential on the upper half-plane, hence has the form $2\pi i f(z)\,dz$ for a weight $2$ modular form $f$, and holomorphy at the cusps forces $f$ to vanish at every cusp. Therefore $\Phi$ is the canonical identification.
[/step]