Let $j$ be Klein's $j$-invariant and let $M$ be a the set of integral $2\times 2$ matrices with determinant $n$, $A$ and $-A$ identified. Why has the polynomial $\Phi(X,Y)\in\mathbb C[X,Y]$ with $$\Phi(X,j(z)) = \prod_{A\in \Gamma\setminus M} (X-j(Az))$$ integral coefficients?
I showed already that there is $\Phi(X,Y)\in\mathbb C[X,Y]$ that satisfies the equation. I did this by regarding the right side as polynomial in $X$ whose coefficients are holomorphic functions in $z$. Then you can see that those coefficients are modular. This implies that they are polynomials in $j$ (each modular function with no poles on $\mathbb H$ is a polynomial in $j$). Is that part correct?