<s>[INST] Proof that $a\nabla^2 u = bu$ is the only homogenous second order 2D PDE unchanged/invariant by rotation Looking for feedback and maybe simpler intuition for my proof of the theorem, shown below
The statement of the theorem:

Theorem 
Among all second-order homogeneous PDEs in two dimensions with constant coefficients, show that the only ones that do not change under a rotation of the coordinate
  system (i.e., are rotationally invariant), have the form
  $$a\nabla^2u = bu $$

Proof:
The general PDE of those conditions is written as:
$$a_1u_{xx} + 2a_2u_{xy} + a_3 u_{yy} + b_1u_x +b_2u_y +cu = 0$$
A counter-clockwise rotation of a point $x,y$ can be given by the rotation matrix is given by the figure below,

where with some basic geometry, we can derive
$$x' = \|{\mathbf{v}}\|\cos\left(\theta + \tan^{-1}\left(\frac{y}{x}\right)\right) = x\cos\theta - y\sin\theta$$
$$y' = \|{\mathbf{v}}\|\sin\left(\theta + \tan^{-1}\left(\frac{y}{x}\right)\right) = x\sin\theta + y\cos\theta$$
This can be summarized in a matrix transformation
$$\begin{bmatrix} x' \\ y' \end{bmatrix} =  
\begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}$$
giving the map:
\begin{align*}
x\mapsto x'=x\cos\theta - y\sin\theta \\
y\mapsto y'=x\sin\theta+y\cos\theta
\end{align*}
From here we find derivatives of our new coordinates:
$$ \frac{\partial x'}{\partial x} = \cos\theta \quad \frac{\partial y'}{\partial x}=\sin\theta$$
$$\frac{\partial x'}{\partial y} = -\sin\theta \quad \frac{\partial y'}{\partial y} = \cos\theta$$
Now the first derivatives of $u(x',y')$ with respect to $x,y$:
\begin{align*}
   u_x = \frac{\partial u}{\partial x} = \frac{\partial u}{\partial x'} \frac{\partial x'}{\partial x} + \frac{\partial u}{\partial y'} \frac{\partial y'}{\partial x} = u_{x'}\cos\theta + u_{y'} \sin\theta \\
   u_y = \frac{\partial u}{\partial x} = \frac{\partial u}{\partial x'} \frac{\partial x'}{\partial y} + \frac{\partial u}{\partial y'} \frac{\partial y'}{\partial y} =  -u_{x'}\sin\theta + u_{y'} \cos\theta
  \end{align*}
And then the second derivatives:
\begin{align*}
   &u_{xx} =  u_{x'x'} \cos^2\theta +2u_{x'y'}\sin\theta\cos\theta+  u_{y'y'} \sin^2\theta \\
   &u_{xy} = -u_{x'x'} \cos\theta\sin\theta - u_{y'x'} \sin^2\theta + u_{x'y'} \cos^2\theta + u_{y'y'}\sin\theta\cos\theta \\
   &u_{yy}=  u_{x'x'} \sin^2\theta -2u_{x'y'}\sin\theta\cos\theta+  u_{y'y'} \cos^2\theta
  \end{align*}
Substituting into the general PDE and rearranging factors of partial derivatives, can be written as
$$ \color{blue}{\widetilde{a_1}u_{x'x'} + \widetilde{a_2}u_{x'y'} + \widetilde{a_3} u_{y'y'} + \widetilde{b_1}u_{x'} +\widetilde{b_2}u_{y'} +\widetilde{c}u = 0} $$
where:
\begin{align*}
    &\widetilde{a_1} =  a_1 \cos^2\theta -2a_2\cos\theta\sin\theta + a_3\sin^2\theta \\
    &\widetilde{a_2} = (a_1-a_3)\sin 2\theta + 2a_2 \cos 2\theta \\
    &\widetilde{a_3} = a_1 \sin^2\theta + 2a_2\sin\theta\cos\theta +a_3\cos^2\theta \\
    &\widetilde{b_1} = b_1\cos\theta - b_2\sin\theta  \\
    &\widetilde{b_2} = b_1\sin\theta + b_2\cos\theta \\
    &\widetilde{c} = c
  \end{align*}
Since we require rotational invariance, the original equation and the transformed PDE must have the same value,namely $0$, on all of $u$. This can be written as
$$\small{a_1u_{xx} + 2a_2u_{xy} + a_3 u_{yy} + b_1u_x +b_2u_y +cu = \widetilde{a_1}u_{x'x'} + \widetilde{a_2}u_{x'y'} + \widetilde{a_3} u_{y'y'} + \widetilde{b_1}u_{x'} +\widetilde{b_2}u_{y'} +\widetilde{c}u} $$
whence:
\begin{align}
   \tag{1} a_1& = a_1 \cos^2\theta -2a_2\cos\theta\sin\theta + a_3\sin^2\theta&\\ 
   \tag{2} 2a_2& = (a_1-a_3)\sin 2\theta + 2a_2 \cos 2\theta& \\
   \tag{3} a_3& = a_1 \sin^2\theta + 2a_2\sin\theta\cos\theta +a_3\cos^2\theta& \\
   \tag{4} b_1& = b_1\cos\theta - b_2\sin\theta&\\
   \tag{5} b_2& = b_1\sin\theta + b_2\cos\theta&\\
   \tag{6} c &= c&
  \end{align}
Excluding the trivial case where $\{a_i\},\{b_i\},c = 0$, we can make several conclusions. Note that the deductions below are made with the understanding that any arbitrary angle $\theta$ must be valid, hence it is erroneous to apply $\theta =0$ in order to reach equality.


*

*can only be true when $a_1=a_3$ and $a_2=0$,

*implies $a_2=0$ and $a_1=a_3$,

*like (1) is only true when  $a_1=a_3$ and $a_2=0$,

*is true when $b_1=b_2=0$

*like (4) is true when $b_1=b_2=0$,

*implies $c\in \mathbb{R}$ is valid.


All together we then know $a1=a3$, $a_2=b_1=b_3=0$ and $c=c$.
Returning to the rotated PDE we now know:
$$\widetilde{a_1} =  a_1,
   \widetilde{a_2} = 0,
   \widetilde{a_3} = a_1 ,
   \widetilde{b_1} = 0,
   \widetilde{b_2} = 0,
   \widetilde{c} = c
  $$
So the PDE under a rotation, $u(x',y')$, becomes
\begin{align*}
   &a_1u_{x'x'} + a_1u_{y'y'} +cu = 0\\
   \Rightarrow& a_1(u_{x'x'}+u_{y'y'}) = -cu \\
   \Rightarrow& a\nabla^2{u} = bu
\end{align*}
where we have chosen $a_1=a, -c = b$ for all $a,b\in\mathbb{R}$. This is in terms of the new rotation $u(x',y')$, and so it remains to show that $\nabla^2{u(x,y)} = \nabla^2{u(x',y')}$, consider from the derivatives before:
        \begin{align*}
   &u_{xx} =  u_{x'x'} \cos^2\theta +2u_{x'y'}\sin\theta\cos\theta+  u_{y'y'} \sin^2\theta \\
   &u_{yy}=  u_{x'x'} \sin^2\theta -2u_{x'y'}\sin\theta\cos\theta+  u_{y'y'} \cos^2\theta
  \end{align*}
Sum them together
$$u_{xx}+u_{yy}=  u_{x'x'} (\sin^2\theta+\cos^2\theta) -2u_{x'y'}\sin\theta\cos\theta+ 2u_{x'y'}\sin\theta\cos\theta+  u_{y'y'} (\cos^2\theta+\sin^2\theta)=u_{x'x'}+u_{y'y'}$$
hence, $\nabla^2{u(x,y)} = \nabla^2{u(x',y')}$ as required.
Therefore a rotation applied to any second order homogeneous 2D PDE with constant coefficents will transform to a PDE of the form $a\nabla^2{u} = bu$ under the rotated coordinate $x',y'$, which we have shown to be equivalent under the regular coordinates $x,y$. This is the only PDE that is invariant under rotation. $$\tag*{$\blacksquare$}$$
Additional Remark
Also was just curious about rotational invariant functions and operators. Anything that solves laplaces equation ($\nabla^2=0$) is called a harmonic function, and satisfies properties such as the mean value property and the maximal principle. I assumed at first harmonic functions meant they were radial, but I think it is more along the lines of being symmetric?
The laplacian is rotationally invariant but the laplace equation has some solutions which are radial (rotationally invariant) and some which are not. I also noticed that the converse is not true, i.e, a radial function does not imply $\nabla^2  = 0$, like $f(x,y)=x^2+y^2$
See bounty remark below
 [/INST] We say that a linear operator $L$ is rotationally invariant if and only if $L$ commutes with the orthogonal group, i.e.
$[L, O] = 0$ for every $O \in \text{O}(n)$. 
Hence what you are proving is that if $L$ is a second order linear operator then 
\begin{align}
LO[f](x) = L[f(O x)] = [Lf](O x) = OL[f](x)
\end{align}
if and only if $L = a\Delta-bI$. Moreover, this is equivalent to showing
\begin{align}
L[f](x, y) = O^{-1}LO[f](x, y)
\end{align}
for every function $f$, that is, $L$ remains fixed under the conjugation action of orthogonal transformations.  
Example:
Let us look at an example. Consider $f(x, y) = x e^y$ and $L=\Delta$. Observe
\begin{align}
O[f]=&\ f(\cos\theta x - \sin\theta y, \sin\theta x+\cos\theta y) \\
=&\ (\cos\theta x-\sin\theta y)e^{\sin\theta x+\cos\theta y}
\end{align}
where
\begin{align}
O
=
\begin{pmatrix}
\cos\theta & -\sin\theta\\
\sin\theta & \cos\theta
\end{pmatrix}.
\end{align}
Then we see that
\begin{align}
g(x, y):=LO[f](x, y)= e^{\sin\theta x+\cos\theta y}(x\cos\theta-y\sin\theta)
\end{align}
and finally
\begin{align}
O^{-1}[g](x, y) =&\ g(\cos\theta x+\sin\theta y, -\sin\theta x+\cos\theta y)\\
=&\ e^{\sin\theta\cos\theta x+\sin^2\theta y-\sin\theta\cos\theta x+\cos^2\theta y}(\cos^2\theta x+\sin\theta \cos\theta y+\sin^2\theta x-\sin\theta\cos\theta y)\\
=&\ xe^y.
\end{align}
Hence
\begin{align}
O^{-1}LO [f](x, y) = xe^y.
\end{align}
Also, note that $\Delta f =x e^y$. Thus, $L[f](x, y) = O^{-1}LO[f](x, y)$. 
Radial Function: In fact the only radial harmonic solution defined on the entire $xy$-plane are constants. This is a simple consequence of the mean-value identity and the maximum principle for harmonic function. Hence $L$ being rotationally invariant doesn't mean
\begin{align}
f(Ox) = f(x) \text{ for all } O \in \text{O}(2)\ \ \implies \ \ \Delta f  =0. 
\end{align} 
Last Remark: Unfortunately, I don't think there are much easier ways to show the only rotationally invariant second order differential operators are given by $L=a\Delta-bI$ other than direct computation. 
</s>