[proofplan]
We prove the identity by computing separately the three factors in the pulled-back Petersson density: the modular-form factor, the imaginary-part factor, and the Lebesgue-area factor. The transformation laws for $f$ and $g$ contribute a factor $|cz+d|^{2k}$. The imaginary part contributes $|cz+d|^{-2k}$, and the hyperbolic area element contributes no remaining factor because the Jacobian $|cz+d|^{-4}$ cancels the square of $\operatorname{Im}(\gamma z)=y|cz+d|^{-2}$.
[/proofplan]
[step:Apply the modular transformation law to the two cusp forms]
Since $f,g \in S_k(SL_2(\mathbb{Z}))$, both are modular forms of weight $k$ for $SL_2(\mathbb{Z})$. Thus, for every $z \in \mathbb{H}$,
\begin{align*}
f(\gamma z) &= (cz+d)^k f(z), \\
g(\gamma z) &= (cz+d)^k g(z).
\end{align*}
Taking the complex conjugate of the second identity gives
\begin{align*}
\overline{g(\gamma z)}
=
\overline{(cz+d)^k g(z)}
=
(\overline{cz+d})^k \overline{g(z)}.
\end{align*}
Multiplying the two identities, we obtain
\begin{align*}
f(\gamma z)\overline{g(\gamma z)}
=
(cz+d)^k(\overline{cz+d})^k f(z)\overline{g(z)}
=
|cz+d|^{2k} f(z)\overline{g(z)}.
\end{align*}
[guided]
The Petersson integrand contains the product $f(\gamma z)\overline{g(\gamma z)}$, so we first compute exactly how this product transforms. Because $f$ and $g$ are cusp forms of weight $k$ for $SL_2(\mathbb{Z})$, they satisfy the same weight-$k$ transformation rule:
\begin{align*}
f(\gamma z) &= (cz+d)^k f(z), \\
g(\gamma z) &= (cz+d)^k g(z).
\end{align*}
The second factor in the Petersson integrand is the complex conjugate of $g(\gamma z)$, so we conjugate the second transformation identity:
\begin{align*}
\overline{g(\gamma z)}
=
\overline{(cz+d)^k g(z)}
=
(\overline{cz+d})^k\overline{g(z)}.
\end{align*}
Multiplying the transformed $f$ factor and the transformed conjugate $g$ factor gives
\begin{align*}
f(\gamma z)\overline{g(\gamma z)}
=
(cz+d)^k(\overline{cz+d})^k f(z)\overline{g(z)}.
\end{align*}
Since $(cz+d)(\overline{cz+d})=|cz+d|^2$, this becomes
\begin{align*}
f(\gamma z)\overline{g(\gamma z)}
=
|cz+d|^{2k}f(z)\overline{g(z)}.
\end{align*}
This is the only place where the modularity of $f$ and $g$ is used.
[/guided]
[/step]
[step:Compute the transformation of the imaginary part]
Let $z=x+iy \in \mathbb{H}$. Since $a,b,c,d \in \mathbb{R}$ and $ad-bc=1$, we compute
\begin{align*}
\gamma z
&=
\frac{az+b}{cz+d}
=
\frac{(az+b)(c\overline{z}+d)}{|cz+d|^2}.
\end{align*}
The imaginary part of the numerator is
\begin{align*}
\operatorname{Im}\bigl((az+b)(c\overline{z}+d)\bigr)
&=
\operatorname{Im}\bigl(ac z\overline{z}+ad z+bc\overline{z}+bd\bigr) \\
&=
ad\,\operatorname{Im}(z)+bc\,\operatorname{Im}(\overline{z}) \\
&=
ad\,y-bc\,y \\
&=
(ad-bc)y \\
&=
y.
\end{align*}
Therefore
\begin{align*}
\operatorname{Im}(\gamma z)
=
\frac{y}{|cz+d|^2},
\end{align*}
and hence
\begin{align*}
\operatorname{Im}(\gamma z)^k
=
y^k |cz+d|^{-2k}.
\end{align*}
[guided]
We next compute the imaginary part because the Petersson density contains the factor $\operatorname{Im}(\gamma z)^k$. Write $z=x+iy$, where $x,y \in \mathbb{R}$ and $y>0$. Since $\gamma \in SL_2(\mathbb{Z})$, its entries are real and satisfy $ad-bc=1$. Rationalising the denominator gives
\begin{align*}
\gamma z
=
\frac{az+b}{cz+d}
=
\frac{(az+b)(c\overline{z}+d)}{|cz+d|^2}.
\end{align*}
The denominator $|cz+d|^2$ is real and positive, so the imaginary part is the imaginary part of the numerator divided by $|cz+d|^2$. Expanding the numerator,
\begin{align*}
(az+b)(c\overline{z}+d)
=
ac z\overline{z}+ad z+bc\overline{z}+bd.
\end{align*}
The terms $ac z\overline{z}=ac|z|^2$ and $bd$ are real. Also $\operatorname{Im}(z)=y$ and $\operatorname{Im}(\overline{z})=-y$. Therefore
\begin{align*}
\operatorname{Im}\bigl((az+b)(c\overline{z}+d)\bigr)
&=
ad\,y+bc(-y) \\
&=
(ad-bc)y \\
&=
y.
\end{align*}
Thus
\begin{align*}
\operatorname{Im}(\gamma z)
=
\frac{y}{|cz+d|^2}.
\end{align*}
Raising both sides to the power $k$ gives
\begin{align*}
\operatorname{Im}(\gamma z)^k
=
y^k |cz+d|^{-2k}.
\end{align*}
This factor will cancel exactly against the $|cz+d|^{2k}$ from the modular transformation law.
[/guided]
[/step]
[step:Compute the Jacobian of the fractional linear map]
Regard $T_\gamma$ as a smooth map
\begin{align*}
T_\gamma : \mathbb{H} &\to \mathbb{H} \\
x+iy &\mapsto x'+iy'.
\end{align*}
As a holomorphic map of the complex variable $z$, its complex derivative is
\begin{align*}
T_\gamma'(z)
=
\frac{a(cz+d)-c(az+b)}{(cz+d)^2}
=
\frac{ad-bc}{(cz+d)^2}
=
\frac{1}{(cz+d)^2}.
\end{align*}
For a holomorphic map, the real Jacobian determinant is the squared modulus of the complex derivative. Hence
\begin{align*}
\left|\det J(T_\gamma)_{(x,y)}\right|
=
|T_\gamma'(z)|^2
=
|cz+d|^{-4}.
\end{align*}
Therefore the Lebesgue area elements satisfy
\begin{align*}
d\mathcal{L}^2(x',y')
=
|cz+d|^{-4}\,d\mathcal{L}^2(x,y).
\end{align*}
[guided]
The area element changes by the real Jacobian determinant of the map $T_\gamma$. We view
\begin{align*}
T_\gamma : \mathbb{H} &\to \mathbb{H} \\
x+iy &\mapsto x'+iy'
\end{align*}
as a real smooth map from the open subset $\mathbb{H}\subset \mathbb{R}^2$ to itself. Since $cz+d\neq 0$ for $z\in \mathbb{H}$, this map is holomorphic on $\mathbb{H}$. Its complex derivative is computed by the quotient rule:
\begin{align*}
T_\gamma'(z)
&=
\frac{a(cz+d)-c(az+b)}{(cz+d)^2} \\
&=
\frac{acz+ad-acz-bc}{(cz+d)^2} \\
&=
\frac{ad-bc}{(cz+d)^2}.
\end{align*}
Because $\gamma \in SL_2(\mathbb{Z})$, we have $ad-bc=1$, so
\begin{align*}
T_\gamma'(z)=\frac{1}{(cz+d)^2}.
\end{align*}
For a holomorphic map $F=u+iv$, the real Jacobian matrix has determinant $|F'(z)|^2$ by the Cauchy-Riemann equations. Applying this to $F=T_\gamma$ gives
\begin{align*}
\left|\det J(T_\gamma)_{(x,y)}\right|
=
|T_\gamma'(z)|^2
=
\left|\frac{1}{(cz+d)^2}\right|^2
=
|cz+d|^{-4}.
\end{align*}
Thus the two-dimensional Lebesgue measure transforms as
\begin{align*}
d\mathcal{L}^2(x',y')
=
|cz+d|^{-4}\,d\mathcal{L}^2(x,y).
\end{align*}
This is the area-scaling factor for the change of variables $z\mapsto \gamma z$.
[/guided]
[/step]
[step:Show that the hyperbolic area factor is invariant]
Using the imaginary-part formula and the Jacobian formula, we obtain
\begin{align*}
\frac{d\mathcal{L}^2(x',y')}{\operatorname{Im}(\gamma z)^2}
&=
\frac{|cz+d|^{-4}\,d\mathcal{L}^2(x,y)}
{\left(y|cz+d|^{-2}\right)^2} \\
&=
\frac{|cz+d|^{-4}\,d\mathcal{L}^2(x,y)}
{y^2|cz+d|^{-4}} \\
&=
\frac{d\mathcal{L}^2(x,y)}{y^2}.
\end{align*}
Thus the hyperbolic area density is invariant under $T_\gamma$.
[guided]
The Petersson integrand contains the hyperbolic area factor
\begin{align*}
\frac{d\mathcal{L}^2(x',y')}{\operatorname{Im}(\gamma z)^2}.
\end{align*}
We now substitute the two formulas already proved:
\begin{align*}
d\mathcal{L}^2(x',y')
=
|cz+d|^{-4}\,d\mathcal{L}^2(x,y)
\end{align*}
and
\begin{align*}
\operatorname{Im}(\gamma z)
=
y|cz+d|^{-2}.
\end{align*}
Therefore
\begin{align*}
\frac{d\mathcal{L}^2(x',y')}{\operatorname{Im}(\gamma z)^2}
&=
\frac{|cz+d|^{-4}\,d\mathcal{L}^2(x,y)}
{\left(y|cz+d|^{-2}\right)^2} \\
&=
\frac{|cz+d|^{-4}\,d\mathcal{L}^2(x,y)}
{y^2|cz+d|^{-4}} \\
&=
\frac{d\mathcal{L}^2(x,y)}{y^2}.
\end{align*}
The factor $|cz+d|^{-4}$ from the Jacobian is exactly cancelled by the square of the imaginary-part scaling. This proves the invariance of the hyperbolic area density.
[/guided]
[/step]
[step:Multiply the transformed factors and cancel the automorphy factors]
Combining the modular-form product, the imaginary-part transformation, and the hyperbolic area invariance, we get
\begin{align*}
&f(\gamma z)\overline{g(\gamma z)}
\operatorname{Im}(\gamma z)^k
\frac{d\mathcal{L}^2(x',y')}{\operatorname{Im}(\gamma z)^2} \\
&\qquad =
\left(|cz+d|^{2k}f(z)\overline{g(z)}\right)
\left(y^k|cz+d|^{-2k}\right)
\frac{d\mathcal{L}^2(x,y)}{y^2} \\
&\qquad =
f(z)\overline{g(z)}y^k
\frac{d\mathcal{L}^2(x,y)}{y^2}.
\end{align*}
This is precisely the claimed invariance of the Petersson inner product integrand.
[/step]
