The Laplacian operator is invariant to $SL_2(\mathbb{R})$

I am reading Iwaniec's book on the spectral analysis of automorphic forms, where I bumped into the following statement in p.20 section 1.6.

Given a function $f:\mathbb{H}\longrightarrow \mathbb{C}$, having continuous second derivatives, and some $g\in SL_2(\mathbb{R})$, then $\Delta(f(gz)) = (\Delta f)(gz)$.

Note that the Laplacian operator we use on the hyperbolic plane $\mathbb{H}$, is $$\Delta = y^2(\dfrac{\partial^2}{\partial x^2} + \dfrac{\partial^2}{\partial y^2})$$

My idea is that in order to prove such a statement, it is enough to prove it over the generators of $SL_2(\mathbb{R})$, i.e., it is enough to check that the identity holds for $g_t(z) = z + t$ for $t\in\mathbb{R}$, and $g^*(z) = -\dfrac{1}{z}$. I found it easy proving the statement for $g=g_t$, however for $g^*$ I didn't understand why it is true that:

$$\Delta(f(-\dfrac{1}{z})) = (\Delta f)(-\dfrac{1}{z})$$

Writing down the formulas for the Laplacian, we get that the LHS is equal to: $$y^2(\dfrac{\partial^2 (f\circ\dfrac{-1}{z})}{\partial x^2} + \dfrac{\partial^2 (f\circ\dfrac{-1}{z})}{\partial y^2})$$ Whereas the RHS equals: $$\Im(-\dfrac{1}{z})^2(\dfrac{\partial^2 f}{\partial x^2}(\dfrac{-1}{z}) + \dfrac{\partial^2 f}{\partial y^2}(\dfrac{-1}{z}))$$

I could continue with the computation, however, simply analyzing the coefficient of $$\dfrac{\partial^2 f}{\partial x^2}(\dfrac{-1}{z})$$ on both sides seems to show me that I am on the wrong path. For the LHS:

$$y^2\cdot\dfrac{\partial^2 f}{\partial x^2}(-\dfrac{1}{z})\dfrac{\partial^2 (-\dfrac{1}{z})}{\partial x^2}$$

For the RHS I obtain:

$$\Im(-\dfrac{1}{z})^2\dfrac{\partial^2 f}{\partial x^2}(\dfrac{-1}{z})$$

Hence the identity could potentially be true if one has the equality: $$\Im(-\dfrac{1}{z})^2 = y^2\dfrac{\partial^2 (-\dfrac{1}{z})}{\partial x^2}$$

However, the RHS is $\dfrac{-2y^2}{z^3}$, whereas the LHS is: $\dfrac{y^2}{(x^2+y^2)^2}$, we cannot expect equality in the general case, so I am confused.

$$(\Delta f)(z)= y^2(f_{xx}(z)+f_{yy}(z))$$

$u(z)=-1/z,g(z)=f(u(z))$

$$g_x(z) = \Re(u_x(z))f_x(u(z))+\Im(u_x(z))f_y(u(z))$$ $$g_{xx}(z) = \Re(u_{xx}(z))f_x(u(z)) + \Re(u_x(z))^2 f_{xx}(u(z))+\Re(u_x(z))\Im(u_x(z)) f_{xy}(u(z))\\ + \Im(u_{xx}(z))f_y(u(z)) + \Re(u_x(z))\Im(u_x(z)) f_{yx}(u(z))+\Im(u_x(z))^2 f_{yy}(u(z)) $$

$$g_y(z) = \Re(u_y(z))f_x(u(z))+\Im(u_y(z))f_y(u(z))$$

$$g_{yy}(z) = \Re(u_{yy}(z))f_x(u(z)) + \Re(u_y(z))^2 f_{xx}(u(z))+\Re(u_y(z))\Im(u_y(z)) f_{xy}(u(z))\\ + \Im(u_{yy}(z))f_y(u(z)) + \Re(u_y(z))\Im(u_y(z)) f_{yx}(u(z))+\Im(u_y(z))^2 f_{yy}(u(z)) $$

The point is that $u_{xx}(z)=-u_{yy}(z)$ and $ \Re(u_x(z))\Im(u_x(z)) = \Re(-i u_y(z))\Im(-iu_y(z))= -\Re(u_y(z))\Im(u_y(z)) $ so you get

$$\boxed{g_{xx}(z)+g_{yy}(z)=(\Re(u_x(z))^2+\Re(u_y(z))^2) f_{xx}(u(z))+ (\Im(u_x(z))^2+\Im(u_y(z))^2) f_{yy}(u(z))}$$ You'll find that $$\Re(u_x(z))^2+\Re(u_y(z))^2=\frac1{|z|^4} = \frac{\Im(-1/z)^2}{y^2}=\Im(u_x(z))^2+\Im(u_y(z))^2$$ and hence $$(\Delta g)(z) = (\Delta f)(-1/z)$$

I think not so many people have ever followed that computation through to the end! :) Ok, so you can see @reuns' good version of it.

Still, I would seriously suggest that we'd want to "somehow" know a-priori that that expression was $SL_2(\mathbb R)$ invariant, rather than be told or have an epiphany that it must be so.

Yes, if we believe already in the existence of an $SL_2(\mathbb R)$-invariant second-order operator (annihilating constants), then already the subgroup $\pmatrix{* & * \cr 0 & *}$ happens to determine the form of it.

But, perhaps, after a certain point in one's life, a more intrinsic and coordinate-independent "characterization" of Casimir/Laplacian might be desirable... so that one can determine what it is in any particular coordinate system, while knowing a-priori that it is invariant. That is, instead of guessing a thing in coordinates and checking its invariance, we characterize an invariant thing and then determine what it is in some given coordinates.

The standard one-line description of a $G$-invariant second-order diffop (Casimir) on a semi-simple real Lie group is obtained from a diagram of $G$-equivariant maps: $$ \mathrm{End}_{\mathbb C}(\mathfrak g) \;\approx\; \mathfrak g\otimes \mathfrak g^* \;\approx\; \mathfrak g\otimes \mathfrak g \;\subset\;A\mathfrak g \;\to\; U\mathfrak g $$ where the first is the general linear algebra isomorphism, the second is the isomorphism of $\mathfrak g$ with its dual via the Killing form, $A\mathfrak g$ is the universal associative algebra generated by $\mathfrak g$ (also called "the tensor algebra"), and the last map is the quotient to the universal enveloping algebra. The identity endomorphism of $\mathfrak g$ is $G$-invariant/equivariant, so its image in $U\mathfrak g$ is, also. That image is the Casimir operator, which descends to the invariant Laplacian on $G/K$...

The identity map in the endomorphism algebra is immediately identifiable with $\sum_i x_i\otimes x_i^*$ for any basis $x_i$ of $\mathfrak g$, with dual basis $x_i^*$. An astute choice of these gives the in-coordinates expression of the invariant Laplacian on $\mathfrak H$... knowing a-priori that whatever expression comes out is $G$-invariant. :)