[proofplan]
We compare the vanishing order of $f$ at $z_0$ with the vanishing order of the pulled-back function $f \circ \gamma$ at $z_0$. The modular transformation law writes $f \circ \gamma$ as $f$ multiplied by the holomorphic nonvanishing factor $(cz+d)^k$, so this multiplication does not change the order at $z_0$. It remains only to check that composing with the local biholomorphism $\gamma$ converts the order of $f \circ \gamma$ at $z_0$ into the order of $f$ at $\gamma z_0$.
[/proofplan]
[step:Reduce the statement to the nonzero modular form case]
If $f$ is the identically zero modular form, then by convention
\begin{align*}
\operatorname{ord}_{\gamma z_0}(f)=\infty=\operatorname{ord}_{z_0}(f),
\end{align*}
so the conclusion follows. Hence assume for the rest of the proof that $f$ is not identically zero.
[/step]
[step:Write the modular transformation law as multiplication by a nonvanishing holomorphic factor]
Let
\begin{align*}
j_\gamma: \mathbb{H} &\to \mathbb{C}\setminus\{0\} \\
z &\mapsto cz+d
\end{align*}
denote the automorphy factor associated to $\gamma$. Since $z_0 \in \mathbb{H}$, the number $cz_0+d$ is nonzero: if $c \neq 0$, then $cz_0+d=0$ would imply $z_0=-d/c \in \mathbb{R}$, contradicting $z_0 \in \mathbb{H}$; if $c=0$, then $ad=1$, so $d=\pm 1$.
Define
\begin{align*}
h: \mathbb{H} &\to \mathbb{C}\setminus\{0\} \\
z &\mapsto j_\gamma(z)^k=(cz+d)^k.
\end{align*}
The modularity of $f \in M_k$ gives, for every $z \in \mathbb{H}$,
\begin{align*}
f(\gamma z)=h(z)f(z).
\end{align*}
Since $h$ is holomorphic and $h(z_0)\neq 0$, multiplication by $h$ does not change the vanishing order at $z_0$. Therefore
\begin{align*}
\operatorname{ord}_{z_0}(f\circ \gamma)=\operatorname{ord}_{z_0}(f).
\end{align*}
[guided]
Let
\begin{align*}
j_\gamma: \mathbb{H} &\to \mathbb{C}\setminus\{0\} \\
z &\mapsto cz+d
\end{align*}
be the automorphy factor. We first verify that this factor does not vanish at the point under consideration. If $c \neq 0$ and $cz_0+d=0$, then $z_0=-d/c$, which is real because $c,d \in \mathbb{Z}$. This contradicts $z_0 \in \mathbb{H}$. If $c=0$, then the condition $\det \gamma=ad-bc=1$ gives $ad=1$, so $d=\pm 1$ and again $cz_0+d=d\neq 0$.
Now define the holomorphic nonvanishing factor
\begin{align*}
h: \mathbb{H} &\to \mathbb{C}\setminus\{0\} \\
z &\mapsto j_\gamma(z)^k=(cz+d)^k.
\end{align*}
The defining transformation law for a weight-$k$ modular form says that
\begin{align*}
f(\gamma z)=h(z)f(z)
\end{align*}
for every $z \in \mathbb{H}$. Near $z_0$, the function $h$ has a convergent [power series](/page/Power%20Series) with nonzero constant term $h(z_0)$. Thus if $f$ has local expansion
\begin{align*}
f(z)=(z-z_0)^m u(z),
\end{align*}
where $m=\operatorname{ord}_{z_0}(f)$ and $u$ is holomorphic with $u(z_0)\neq 0$, then
\begin{align*}
h(z)f(z)=(z-z_0)^m h(z)u(z),
\end{align*}
and $h(z_0)u(z_0)\neq 0$. Hence $h f$ has the same vanishing order as $f$ at $z_0$. Since $f\circ\gamma=h f$, this proves
\begin{align*}
\operatorname{ord}_{z_0}(f\circ \gamma)=\operatorname{ord}_{z_0}(f).
\end{align*}
[/guided]
[/step]
[step:Identify the order of the pullback with the order at the image point]
Define
\begin{align*}
\varphi_\gamma: \mathbb{H} &\to \mathbb{H} \\
z &\mapsto \gamma z=\frac{az+b}{cz+d}.
\end{align*}
This map is biholomorphic, with inverse $\varphi_{\gamma^{-1}}$. Its complex derivative at $z_0$ is
\begin{align*}
\varphi_\gamma'(z_0)=\frac{ad-bc}{(cz_0+d)^2}=\frac{1}{(cz_0+d)^2}\neq 0.
\end{align*}
Let $w_0=\gamma z_0$. Since $f$ is holomorphic and not identically zero, there are an integer $m\geq 0$, a neighbourhood $V \subset \mathbb{H}$ of $w_0$, and a [holomorphic function](/page/Holomorphic%20Function)
\begin{align*}
u: V &\to \mathbb{C}
\end{align*}
such that $u(w_0)\neq 0$ and
\begin{align*}
f(w)=(w-w_0)^m u(w)
\end{align*}
for every $w \in V$. By definition, $m=\operatorname{ord}_{w_0}(f)=\operatorname{ord}_{\gamma z_0}(f)$.
For $z$ sufficiently close to $z_0$, we have $\varphi_\gamma(z)\in V$, and hence
\begin{align*}
(f\circ\gamma)(z)
&=f(\varphi_\gamma(z)) \\
&=(\varphi_\gamma(z)-w_0)^m u(\varphi_\gamma(z)).
\end{align*}
Because $\varphi_\gamma'(z_0)\neq 0$, the holomorphic function
\begin{align*}
q: U &\to \mathbb{C} \\
z &\mapsto
\begin{cases}
\dfrac{\varphi_\gamma(z)-\varphi_\gamma(z_0)}{z-z_0}, & z\neq z_0,\\
\varphi_\gamma'(z_0), & z=z_0
\end{cases}
\end{align*}
is defined on a sufficiently small neighbourhood $U \subset \mathbb{H}$ of $z_0$ and satisfies $q(z_0)\neq 0$. Thus
\begin{align*}
\varphi_\gamma(z)-w_0=(z-z_0)q(z),
\end{align*}
and therefore
\begin{align*}
(f\circ\gamma)(z)=(z-z_0)^m q(z)^m u(\varphi_\gamma(z)).
\end{align*}
The factor $q(z)^m u(\varphi_\gamma(z))$ is holomorphic near $z_0$ and has nonzero value at $z_0$. Hence
\begin{align*}
\operatorname{ord}_{z_0}(f\circ\gamma)=m=\operatorname{ord}_{\gamma z_0}(f).
\end{align*}
[guided]
We now justify the coordinate-change part of the argument. Define
\begin{align*}
\varphi_\gamma: \mathbb{H} &\to \mathbb{H} \\
z &\mapsto \gamma z=\frac{az+b}{cz+d}.
\end{align*}
This is the fractional linear transformation associated to $\gamma$. Its inverse is the fractional linear transformation associated to $\gamma^{-1}$, so $\varphi_\gamma$ is biholomorphic. At the point $z_0$, its derivative is
\begin{align*}
\varphi_\gamma'(z_0)=\frac{ad-bc}{(cz_0+d)^2}=\frac{1}{(cz_0+d)^2}.
\end{align*}
The denominator is nonzero by the previous step, so $\varphi_\gamma'(z_0)\neq 0$.
Set $w_0=\gamma z_0$. By the definition of vanishing order for a nonzero holomorphic function, there are an integer $m\geq 0$, a neighbourhood $V \subset \mathbb{H}$ of $w_0$, and a holomorphic function
\begin{align*}
u: V &\to \mathbb{C}
\end{align*}
with $u(w_0)\neq 0$ such that
\begin{align*}
f(w)=(w-w_0)^m u(w)
\end{align*}
for every $w \in V$. This integer is exactly
\begin{align*}
m=\operatorname{ord}_{w_0}(f)=\operatorname{ord}_{\gamma z_0}(f).
\end{align*}
For $z$ in a sufficiently small neighbourhood $U \subset \mathbb{H}$ of $z_0$, the point $\varphi_\gamma(z)$ lies in $V$. Substituting $w=\varphi_\gamma(z)$ into the local factorization of $f$ gives
\begin{align*}
(f\circ\gamma)(z)
&=f(\varphi_\gamma(z)) \\
&=(\varphi_\gamma(z)-w_0)^m u(\varphi_\gamma(z)).
\end{align*}
The remaining question is whether $\varphi_\gamma(z)-w_0$ vanishes to first order in $z-z_0$. Since $\varphi_\gamma'(z_0)\neq 0$, it does. Explicitly, define
\begin{align*}
q: U &\to \mathbb{C} \\
z &\mapsto
\begin{cases}
\dfrac{\varphi_\gamma(z)-\varphi_\gamma(z_0)}{z-z_0}, & z\neq z_0,\\
\varphi_\gamma'(z_0), & z=z_0.
\end{cases}
\end{align*}
The removable singularity at $z_0$ is filled by the derivative value, so $q$ is holomorphic on $U$, and $q(z_0)=\varphi_\gamma'(z_0)\neq 0$. Therefore
\begin{align*}
\varphi_\gamma(z)-w_0=(z-z_0)q(z).
\end{align*}
Substituting this identity into the expression for $f\circ\gamma$ yields
\begin{align*}
(f\circ\gamma)(z)=(z-z_0)^m q(z)^m u(\varphi_\gamma(z)).
\end{align*}
The factor $q(z)^m u(\varphi_\gamma(z))$ is holomorphic near $z_0$, and its value at $z_0$ is
\begin{align*}
q(z_0)^m u(w_0),
\end{align*}
which is nonzero. Hence $f\circ\gamma$ vanishes to order $m$ at $z_0$, so
\begin{align*}
\operatorname{ord}_{z_0}(f\circ\gamma)=m=\operatorname{ord}_{\gamma z_0}(f).
\end{align*}
[/guided]
[/step]
[step:Combine the two order identities]
The previous two steps give
\begin{align*}
\operatorname{ord}_{\gamma z_0}(f)
=
\operatorname{ord}_{z_0}(f\circ\gamma)
=
\operatorname{ord}_{z_0}(f).
\end{align*}
This is the desired invariance of vanishing order along the $SL_2(\mathbb{Z})$-orbit of $z_0$.
[/step]
