[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.[/step]

[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: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*}[/step]

[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: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))$.[/step]

[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: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.[/step]

[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]