[proofplan]
We compute the Euler characteristic of $X_0(N)$ by viewing $\Gamma_0(N)\backslash \mathbb{H}$ as a hyperbolic orbifold and then adding its cusps. The quotient by $SL_2(\mathbb{Z})$ has orbifold signature $(0;2,3;1)$, so its orbifold Euler characteristic is $-1/6$. Since $\Gamma_0(N)$ has index $\mu_N$, the orbifold Euler characteristic scales by $\mu_N$. Finally, rewriting the orbifold Euler characteristic of $X_0(N)$ in terms of its genus, elliptic points, and cusps gives the stated formula.
[/proofplan]
[step:Identify the relevant orbifold cover]
Let $\mathbb{H} := \{z \in \mathbb{C} : \operatorname{Im}(z) > 0\}$ be the upper half-plane. The group $SL_2(\mathbb{Z})$ acts on $\mathbb{H}$ by fractional linear transformations,
\begin{align*}
\begin{pmatrix}
a & b\\
c & d
\end{pmatrix}
\cdot z
=
\frac{az+b}{cz+d}.
\end{align*}
Define the projective modular group
\begin{align*}
PSL_2(\mathbb{Z}) := SL_2(\mathbb{Z}) / \{\pm I\},
\end{align*}
where $I$ is the $2 \times 2$ identity matrix. Since $-I \in \Gamma_0(N)$, the effective acting group is the image of $\Gamma_0(N)$ under the quotient homomorphism $SL_2(\mathbb{Z}) \to PSL_2(\mathbb{Z})$. The index of this image in $PSL_2(\mathbb{Z})$ is still
\begin{align*}
\mu_N = [SL_2(\mathbb{Z}) : \Gamma_0(N)],
\end{align*}
because quotienting both groups by the common central subgroup $\{\pm I\}$ preserves the coset set.
Thus the natural map of compactified modular curves
\begin{align*}
\pi: X_0(N) \longrightarrow X(1)
\end{align*}
has degree $\mu_N$, where $X(1)$ is the compactification of $SL_2(\mathbb{Z})\backslash \mathbb{H}$.
[guided]
The matrix $-I$ fixes every point of $\mathbb{H}$ under the fractional linear action. Therefore the effective group of transformations is the projective modular group
\begin{align*}
PSL_2(\mathbb{Z}) := SL_2(\mathbb{Z})/\{\pm I\}.
\end{align*}
Because $-I$ belongs to $\Gamma_0(N)$ for every $N$, the quotient homomorphism $SL_2(\mathbb{Z}) \to PSL_2(\mathbb{Z})$ sends $\Gamma_0(N)$ to a subgroup whose cosets are naturally the same as the cosets of $\Gamma_0(N)$ in $SL_2(\mathbb{Z})$. Hence passing from $SL_2(\mathbb{Z})$ to $PSL_2(\mathbb{Z})$ does not change the subgroup index, and the degree of the covering of modular orbifolds
\begin{align*}
\pi: X_0(N) \longrightarrow X(1)
\end{align*}
is precisely
\begin{align*}
\mu_N = [SL_2(\mathbb{Z}) : \Gamma_0(N)].
\end{align*}
This is the point where the index $\mu_N$ enters the genus formula.
[/guided]
[/step]
[step:Compute the orbifold Euler characteristic of the base curve]
The compact modular curve $X(1)$ has underlying Riemann surface of genus $0$, one elliptic point of order $2$, one elliptic point of order $3$, and one cusp. For any finite-area modular orbifold $Y$, let $\chi_{\mathrm{orb}}(Y)$ denote its orbifold Euler characteristic, computed from the underlying compact genus, the elliptic stabilizer orders, and the number of cusps. Therefore
\begin{align*}
\chi_{\mathrm{orb}}(X(1))
&=
2 - 2\cdot 0
-\left(1-\frac{1}{2}\right)
-\left(1-\frac{1}{3}\right)
-1 \\
&=
2-\frac{1}{2}-\frac{2}{3}-1 \\
&=
-\frac{1}{6}.
\end{align*}
Here we use the standard orbifold Euler characteristic formula for finite-area modular orbifolds with cusps.
Since orbifold Euler characteristic multiplies under finite orbifold covers, and since the previous step identified $X_0(N)\to X(1)$ as a finite orbifold cover of degree $\mu_N$, this gives
\begin{align*}
\chi_{\mathrm{orb}}(X_0(N))
=
\mu_N \chi_{\mathrm{orb}}(X(1))
=
-\frac{\mu_N}{6}.
\end{align*}
[guided]
The quotient $SL_2(\mathbb{Z})\backslash \mathbb{H}$ has one elliptic orbit of stabilizer order $2$, one elliptic orbit of stabilizer order $3$, and one cusp. After compactification, the underlying compact Riemann surface is $X(1)\cong \mathbb{P}^1(\mathbb{C})$, which has genus $0$.
For a finite-area modular orbifold $Y$ with underlying compact genus $g$, elliptic points of orders $m_1,\dots,m_r$, and $c$ cusps, define $\chi_{\mathrm{orb}}(Y)$ to be its orbifold Euler characteristic. The standard formula is
\begin{align*}
\chi_{\mathrm{orb}}(Y)
=
2-2g-\sum_{j=1}^{r}\left(1-\frac{1}{m_j}\right)-c.
\end{align*}
Applying this with $g=0$, elliptic orders $2$ and $3$, and $c=1$, we get
\begin{align*}
\chi_{\mathrm{orb}}(X(1))
&=
2-\left(1-\frac{1}{2}\right)-\left(1-\frac{1}{3}\right)-1 \\
&=
2-\frac{1}{2}-\frac{2}{3}-1 \\
&=
-\frac{1}{6}.
\end{align*}
The previous step identified the map $X_0(N)\to X(1)$ as a finite orbifold cover of degree $\mu_N$. Orbifold Euler characteristic is multiplicative under finite orbifold covers, so the required hypothesis for multiplicativity is satisfied, and
\begin{align*}
\chi_{\mathrm{orb}}(X_0(N))
=
\mu_N \chi_{\mathrm{orb}}(X(1))
=
-\frac{\mu_N}{6}.
\end{align*}
[/guided]
[/step]
[step:Express the same Euler characteristic using the signature of $X_0(N)$]
The modular curve $X_0(N)$ has underlying compact Riemann surface of genus $g(X_0(N))$, has $e_2(N)$ elliptic points of order $2$, has $e_3(N)$ elliptic points of order $3$, and has $c(N)$ cusps. Therefore
\begin{align*}
\chi_{\mathrm{orb}}(X_0(N))
&=
2-2g(X_0(N))
-
e_2(N)\left(1-\frac{1}{2}\right)
-
e_3(N)\left(1-\frac{1}{3}\right)
-
c(N) \\
&=
2-2g(X_0(N))
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N).
\end{align*}
Combining this with $\chi_{\mathrm{orb}}(X_0(N))=-\mu_N/6$ gives
\begin{align*}
2-2g(X_0(N))
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N)
=
-\frac{\mu_N}{6}.
\end{align*}
[guided]
Now we compute the same invariant from the geometry of $X_0(N)$ itself. By definition, $e_2(N)$ counts the elliptic points whose stabilizer has order $2$, $e_3(N)$ counts the elliptic points whose stabilizer has order $3$, and $c(N)$ counts the cusps.
The orbifold Euler characteristic formula therefore gives
\begin{align*}
\chi_{\mathrm{orb}}(X_0(N))
&=
2-2g(X_0(N))
-
e_2(N)\left(1-\frac{1}{2}\right)
-
e_3(N)\left(1-\frac{1}{3}\right)
-
c(N) \\
&=
2-2g(X_0(N))
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N).
\end{align*}
But the previous step computed the same quantity from the covering $X_0(N)\to X(1)$:
\begin{align*}
\chi_{\mathrm{orb}}(X_0(N))
=
-\frac{\mu_N}{6}.
\end{align*}
Equating the two expressions yields
\begin{align*}
2-2g(X_0(N))
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N)
=
-\frac{\mu_N}{6}.
\end{align*}
[/guided]
[/step]
[step:Solve the Euler characteristic identity for the genus]
Starting from
\begin{align*}
2-2g(X_0(N))
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N)
=
-\frac{\mu_N}{6},
\end{align*}
we move all terms except the genus term to the right-hand side:
\begin{align*}
2g(X_0(N))
=
2+\frac{\mu_N}{6}
-\frac{e_2(N)}{2}
-\frac{2e_3(N)}{3}
-c(N).
\end{align*}
Dividing by $2$ gives
\begin{align*}
g(X_0(N))
=
1+\frac{\mu_N}{12}
-\frac{e_2(N)}{4}
-\frac{e_3(N)}{3}
-\frac{c(N)}{2}.
\end{align*}
This is the desired formula.
[/step]
