Show Laplace operator is rotationally invariant

Here is an answer along the lines of @TeeJay, but with a few more details and in general $\mathbb R^n$ space. So, we need to show that if $u$ is harmonic ($\Delta_y u(y)=0$) then $v(x)=u(Mx)$ is also harmonic ($\Delta_x v(x)=0$) for an orthogonal $M$ (i.e., $M^\top=M^{-1}$).

First I note that if $D^2_y u(y)$ denotes the Hessian matrix, then $$ \Delta_y u(y)={\rm tr}\, D^2_y u(y). $$ Note also that $D^2=\nabla\nabla^T$, where $\nabla$ is the column vector of partial derivatives.

Now the only thing we need is the fact that $\nabla_x u(Mx)=M^\top \nabla_y u(y)|_{y=Mx}$ and that ${\rm tr}\,(AB)={\rm tr}\,(BA)$.

$$ \Delta_x v(x)={\rm tr}(D^2_x v(x))={\rm tr}(D^2_x u(Mx))=\\ {\rm tr}(\nabla_x\nabla_x^Tu(Mx))={\rm tr}(\nabla_x(\nabla_x u(Mx))^\top)=\\ {\rm tr}(\nabla_x(M^\top \nabla_y u(y))^\top)={\rm tr}(\nabla_x(\nabla_y^\top u(y))M)=\\ {\rm tr}(M^\top\nabla_y\nabla_y^\top u(y)M)={\rm tr}(\nabla_y\nabla_y^\top u(y))\\ ={\rm tr}(D^2_y u(y))=\Delta_y u(y)=0 $$ as required.

$$ \frac{\partial f}{\partial y} = \frac{\partial f}{\partial u} \frac{\partial u}{\partial y} + \frac{\partial f}{\partial v} \frac{\partial v}{\partial y} = \frac{\partial f}{\partial u} \sin \theta+ \frac{\partial f}{\partial v} \cos \theta$$

$$\frac{\partial^2f}{\partial y^2}=\frac{\partial}{\partial y} \left ( \frac{\partial f}{\partial y} \right)=\frac{\partial}{\partial y} \left ( \frac{\partial f}{\partial u} \sin \theta + \frac{\partial f}{\partial v} \cos \theta \right)= \frac{\partial}{\partial y} \frac{\partial f}{\partial u} \sin \theta + \frac{\partial}{\partial y} \frac{\partial f}{\partial v} \cos \theta$$

Now the problem is to compute $\frac{\partial}{\partial y} \frac{\partial f}{\partial u}$ and $\frac{\partial } {\partial y} \frac{\partial f}{\partial v}$. We can think of $\frac{\partial}{\partial y} \frac{\partial f}{\partial u}$ as $\frac{\partial}{\partial u} \frac{\partial f}{\partial y}$ and similarly for the latter.

$$\frac{\partial}{\partial u} \frac{\partial f}{\partial y} = \frac{\partial}{\partial u} \left( \frac{\partial f}{\partial u} \sin \theta + \frac{\partial f}{\partial v} \cos \theta \right) = \frac{\partial^2 f}{\partial u^2} \sin \theta + \frac{\partial f}{\partial u \partial v} \cos \theta$$

$$\frac{\partial}{\partial v} \frac{\partial f}{\partial y} = \frac{\partial}{\partial v} \left( \frac{\partial f}{\partial u} \sin \theta + \frac{\partial f}{\partial v} \cos \theta \right) = \frac{\partial^2 f}{\partial u \partial v} \sin \theta + \frac{\partial^2 f}{\partial v^2} \cos \theta$$

Now plug this back in to $\frac{\partial^2f}{\partial y^2}$ giving,

$$\frac{\partial^2f}{\partial y^2}= \left (\frac{\partial^2 f}{\partial u^2} \sin \theta + \frac{\partial f}{\partial u \partial v} \cos \theta \right) \sin \theta + \left (\frac{\partial^2 f}{\partial u \partial v} \sin \theta + \frac{\partial^2 f}{\partial v^2} \cos \theta \right) \cos \theta $$

Repeat the process for $\frac{\partial^2 f}{\partial x^2}$, add together, and you will arrive at the desired result.

Off the top of my head, I think an easier way to do this would be to write it all in vector/matrix notation. Under a rotation, $\vec\nabla f\rightarrow R\cdot \vec\nabla f$ with $R$ being your rotation matrix. Then the Laplacian would transform like $\nabla^{2}f=\vec\nabla\cdot \vec\nabla f\rightarrow \vec\nabla\cdot R^{T}R\cdot \vec\nabla f$. Rotation matrices satisfy $R^{T}R=1$ so that should do it. There is a bit more rigor to this but this should be a good starting point.