Let $X_0(N) = \Gamma_0(N) / (\mathbb{H} \cup \mathbb{P}^1(\mathbb{Q}))$. A lecture note (p. 2) lists the following (very easy!) formula for the corresponding cusps:
$$\nu_\infty = \sum_{d\mid N} \varphi(\gcd(d,N/d)).$$
Unfortunately, there is no source or derivation given for this formula. I know what a cusp is and that the set of cusps is finite for all subcongruence groups of $\operatorname{SL}_2(\mathbb{Z})$. However, I do not grasp why this formula for $\nu_\infty$ is correct.
See Section 3.8 of Diamond and Shurman's "A first course in modular forms", page 103.
Alternatively, see Prop. 1.43 in page 24 of Shimura's "Introduction to the arithmetic theory of automorphic functions".
This is proposition 2.2 in Yuri Manin's paper 'Parabolic points and zeta functions of modular curves'. He proves a bit more namely that the map $$ \frac{u}{v\delta}\mapsto [\delta,uv\mod \gcd(\delta,N\delta^{-1})], \quad \infty\mapsto[N,1 ] $$ is a bijection, the domain being $\Gamma_0(N)\backslash P^1(\mathbb{Q})$ the cusps of $\Gamma_0(N)$ and the codomain the pairs $[\delta,x\mod\gcd(\delta,N\delta^{-1})]$, where $\delta$ is a divisor of $N$ and $x$ is invertible modulo $\gcd(\delta,N\delta^{-1}]$. $\delta$ is chosen so that $\gcd(u,v\delta)=\gcd(v,N\delta^{-1})=1$. The second set in this bijection obviously has the right cardinality.
It's not hard to see the map in the second answer is injective, as follows.
It suffices to show that if two element
\begin{equation} \frac{u_1}{v_1\delta}, \frac{u_2}{v_2\delta} \end{equation}
satisfy
\begin{equation} u_1v_1\equiv u_2v_2 (\text{mod gcd}(\delta,\frac N\delta)), \end{equation}
then I can find a matrix in $\Gamma_0(N)$ that moves the first element to the second.(where of course, $\delta|N$, $u_i$ is coprime to $v_i\delta$)
Firstly, WLOG I can assume that $v_1=v_2=1$, and I am handling
\begin{equation} \frac{u_1}{\delta}, \frac{u_2}{\delta}, \end{equation}
because, by B$\acute{e}$zout lemma, for any $\frac{u}{v\delta}$, I can find a matrix $\left(\begin{matrix}a&b\\ kN&d\\\end{matrix}\right)$ such that
\begin{equation} kNu+dv\delta=\delta. \end{equation}
then apply this matrix, I will get
\begin{equation} \left(\begin{matrix}a&b\\ kN&d\\\end{matrix}\right)\frac{u}{v\delta}=\frac{au+bv\delta}{\delta}. \end{equation}
Secondly, I can use matrix $\left(\begin{matrix}1&b\\ 0&1\\\end{matrix}\right)$,
to plus/minus $u_1$ or $u_2$, with any multiple of $\delta$, so all I need to care is the class of $u_1$ and $u_2$ mod $\delta$.
So the only thing I need to do, is to find a matrix $\left(\begin{matrix}a&b\\ kN&d\\\end{matrix}\right)$ in $\Gamma_0(N)$ that moves $\frac{u_1}{\delta}$ to $\frac{u_2+?\delta}{\delta}$, under the condition $u_1\equiv u_2$ (mod $gcd(\delta,\frac{N}{\delta})$).
I need to satisfy these conditions:
(1) \begin{equation} au_1+b\delta\equiv u_2 \text{ (mod} \delta) \end{equation}
(2) \begin{equation} kNu_1+d\delta=\delta \end{equation}
(3) \begin{equation} ad-kNb=1. \end{equation}
This is basically B$\acute{e}$zout lemma, really, it's complicated but elementary.
Let me explain how to quickly see this.
(1) is equivalent to $a\equiv [\frac{u_2}{u_1}]$ (mod $\delta$), where $[\frac{u_2}{u_1}]$ is an integer, and $[\frac{u_2}{u_1}]\equiv 1$ (mod $(\delta,\frac{N}{\delta})$).
(2) is equivalent to $d=1-k\frac{N}{\delta}u_1$.
(3) is equivalent to $kN|ad-1$, so that I can let $b=\frac{ad-1}{kN}$.
Let \begin{equation}a=[\frac{s_2}{s_1}]+l\delta+nN\end{equation} where $l,n$ are two new variables to be determined. (1) is now automatically satisfied.
Let \begin{equation} d=1-k\frac{N}{\delta}u_1, \end{equation} so (2) is also equivalently satisfied.
Now we need to carefully chose $l,k,n$ such that (3) holds:
\begin{equation} ad-1= ([\frac{u_2}{u_1}]+l\delta+nN)(1-k\frac{N}{\delta}u_1)-1 \end{equation}
By B$\acute{e}$zout lemma, I can firstly carefully choose $l$ and $k$ such that $ad-1$ is divisible by $N$, (and now k is fixed!) and then I can carefully choose $n$ such that $ad-1$ is divisible by $kN$.
So now I have completed it.
