There is such a theory, but the analytic object that the forms live on is an analytified modular curve, not simply $\mathbb{C}_p$ (though there is a "$p$-adic upper-half plane" that can be used to uniformize some similar moduli spaces, but as far as I know not the usual modular curves).
Basically, if $f$ is a classical modular form of some weight $k$, $f$ can be realized as a section of a sheaf $\omega^{\otimes k}$ on a complex-analytic modular curve such as $X_1(N)$ obtained via quotient from the complex upper-half plane. These curves and sheaves have algebraic models defined over $\mathbb{Q}$, and (under mild hypotheses) the form $f$ actually arises from a section of the associated sheaf $\omega^{\otimes k}$ on the modular curve defined over $\mathbb{Q}$ (or some finite extension).
Now you can go in another direction and consider the rigid-analytic space over $\mathbb{Q}_p$ associated to the smooth algebraic curve $X_1(N)$ and your form $f$ gives rise to a $p$-adic analytic object on this curve. Now one can play games like considering subspaces obtained by removing disks around supersingular points to obtain general $p$-adic modular forms (such as limits of classical forms of varying weight) and overconvergent modular forms.
For what it's worth, these curves are $p$-adic analytic moduli spaces, and this point of view on $p$-adic modular forms essentially differs from Katz's by thinking about rigid spaces as opposed to formal schemes (the Raynaud point of view on $p$-adic analytic geometry).
Look at the papers of Coleman (such as his $p$-adic Banach spaces paper or the eigencurve paper with Mazur) if you would like to read more on this point of view.
Algebraicity is proven in appendix one of [1]. The function considered there is
$$ \psi(\tau) = \frac{3E_4(\tau)}{2E_6(\tau)} (E_2(\tau) - \frac{3}{\pi \rm{Im} \tau}) = \frac{3}{2} s_2(\tau).$$
The proof is by establishing, for quadratic irrationals $\tau$, the identity
$$ \psi(\tau) = 9j(\tau)\gamma + \frac{3(7j(\tau)-6912)}{2(j(\tau)-1728)}$$
where $\gamma$ is a rational function in the coefficients of the Taylor series of a certain polynomial $\Phi(X,Y)$ around the point $(j(\tau),j(\tau))$ where it vanishes. Then
$$M(\tau) = 1728 s_2(\tau) (j(\tau)-1) = 1728(j(\tau)-1) ( 6j(\tau)\gamma + \frac{(7j(\tau) - 4\cdot1728)}{j(\tau)-1728}.$$
There is also a different identity for calculations, for when $\tau$ is not equivalent to $i$ or $\rho$, and satisfies
$$ A \tau^2 + B \tau + C = 0$$
with $A,B,C$ coprime positive integers:
$$\psi(\tau) = \frac{-g_2 S}{Cg_3 (2A+B\tau)}$$
where $S$ is the sum of $\wp(z)$ as $z$ ranges over $C\tau$-torsion points of $\mathbb{C}/\langle 1,\tau \rangle$.
Masser gives a table of $\psi(\tau)$, for a list of generators of rings of integers of imaginary quadratic fields with class number 1. They are all rational numbers, and one can indeed check that
$$ M(\frac{-1+\sqrt{-163}}{2}) = 223263987730882560,$$
which agrees with the factorization in the question. Also
$$\begin{align}M(\frac{-1+\sqrt{-67}}{2})&=-112852776960=-2^{11}\cdot 3^5 \cdot 5 \cdot 7 \cdot 11 \cdot 19 \cdot 31\\
M(\frac{-1+\sqrt{-43}}{2})&=-627056640=-2^{13}\cdot 3^7 \cdot 5 \cdot 7\\
M(\frac{-1+\sqrt{-27}}{2})&=-7772161=-1181\cdot 6581\\
M(\frac{-1+\sqrt{-19}}{2})&=-497664=-2^{11} \cdot 3^5\\
M(\frac{-1+\sqrt{-11}}{2})&=-14336=-2^{11} \cdot 7\end{align},$$
etc.
Probably one can prove $M(\tau)$ is an algebraic integer by going through Masser's calculations and clearing denominators.
[1] Masser, David W. Elliptic Functions and Transcendence
Best Answer
The constant $\alpha$ in your question can be in fact written explicitly as $(k)_n/12^n$, where $(a)_n=\Gamma(a+n)/\Gamma(n)$ is the Pochhammer symbol (shifted factorial) and $k$ denotes the (even) weight of the corresponding Eisenstein series.
Your observation is indeed related to the Rankin--Cohen brakets; see Section 5.2 in [D. Zagier, Elliptic modular forms and their applications, The 1-2-3 of modular forms, Universitext (Springer, Berlin, 2008), pp. 1–-103]. Preserving the notation $D$ of Zagier's lectures for your differential operator and picking a modular form $f$ of weight $k$, one can show that $D^nf$ transforms under the modular group as $$ D^nf\biggl(\frac{a\tau+b}{c\tau+d}\biggr) =\sum_{r=0}^n\binom{n}{r}\frac{(k+r)_{n-r}}{(2\pi i)^{n-r}} c^{n-r}(c\tau+d)^{k+n+r}D^rf(\tau), $$ by the induction on $n\ge 0$. In addition, the function $E_2(\tau)$ transforms as $$ E_2\biggl(\frac{a\tau+b}{c\tau+d}\biggr) =\frac{12c(c\tau+d)}{2\pi i}+(c\tau+d)^2E_2(\tau). $$ Therefore, it remains to verify that the difference $$ g_n=D^nE_k-\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom{n}{r}E_{k+2n-2r}E_2^r $$ satisfies $$ g_n\biggl(\frac{a\tau+b}{c\tau+d}\biggr)=(c\tau+d)^{k+2n}g_n(\tau). $$
Technicalities. Indeed, I found the remaining details quite boring, but going through my yesterday writing I have realised that your expectation fails already for $D^5E_4$ and $D^4E_6$. Here is my explanation why.
Because $g_n(\tau)$ is a $q$-series, so it is invariant under $\tau\mapsto\tau+1$, we can restrict to verifying the claim under the transformation $\tau\mapsto-1/\tau$ (that is, $a=d=0$, $b=-1$, and $c=1$). Then (we take $s=n-r$ in the above formula) $$ D^nE_k(-1/\tau) =\sum_{s=0}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}D^{n-s}E_k(\tau) $$ and $$ \begin{aligned} & \frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nrE_{k+2n-2r}E_2^r\bigg|_{\tau\mapsto-1/\tau} \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nr\tau^{k+2n-r}E_{k+2n-2r}\biggl(\tau E_2+\frac{12}{2\pi i}\biggr)^r \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nr\tau^{k+2n-r}E_{k+2n-2r}\sum_{s=0}^r\binom rs\frac{12^s}{(2\pi i)^s}\tau^{r-s}E_2^{r-s} \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{s=0}^n\binom ns\tau^{r-s}\frac{12^s}{(2\pi i)^s}\tau^{k+2n-s} \sum_{r=s}^n(-1)^{n-r}\binom{n-s}{r-s}E_{k+2n-2r}E_2^{r-s}. \end{aligned} $$ Subtracting the latter from the former we obtain $$ \begin{aligned} g_n(-1/\tau) &=\sum_{s=0}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}g_{n-s}(\tau) \cr &=\tau^{k+2n}g_n(\tau)+\sum_{s=1}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}g_{n-s}(\tau). \end{aligned} $$ Therefore, $g_n(-1/\tau)=\tau^{k+2n}g_n(\tau)$, hence $g_n(\tau)$ is a modular form (of weight $k+2n$), if and only if the additional sum over $s$ vanishes, that is, $g_{n-s}=0$ for $s=1,\dots,n$. The latter however does not happen when $k+2n>12$.