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
 A: 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. 
A: Perhaps writing the real variables $x$ and $y$ as complex variables $z$ and $\bar{z}$ could provide some information as expected.
Define
\begin{align}
z&=x+iy,\\
\bar{z}&=x-iy,
\end{align}
which yields
\begin{align}
\frac{\partial}{\partial x}&=\frac{\partial}{\partial z}+\frac{\partial}{\partial\bar{z}},\\
\frac{\partial}{\partial y}&=i\left(\frac{\partial}{\partial z}-\frac{\partial}{\partial\bar{z}}\right).
\end{align}
Thanks to these relations, we have
\begin{align}
u_x&=u_z+u_{\bar{z}},\\
u_y&=i\left(u_z-u_{\bar{z}}\right),\\
u_{xx}&=u_{zz}+2u_{z\bar{z}}+u_{\bar{z}\bar{z}},\\
u_{xy}&=i\left(u_{zz}-u_{\bar{z}\bar{z}}\right),\\
u_{yy}&=-\left(u_{zz}-2u_{z\bar{z}}+u_{\bar{z}\bar{z}}\right).
\end{align}
As a consequence,
$$
a_1u_{xx}+2a_2u_{xy}+a_3u_{yy}+b_1u_x+b_2u_y+cu=0
$$
is equivalent to
\begin{equation}
\left(a_1+2ia_2-a_3\right)u_{zz}+2\left(a_1+a_3\right)u_{z\bar{z}}+\left(a_1-2ia_2-a_3\right)u_{\bar{z}\bar{z}}+\left(b_1+ib_2\right)u_z+\left(b_1-ib_2\right)u_{\bar{z}}+cu=0.\tag{1}
\end{equation}
Now, for rotational transformation, we have
$$
z\to e^{i\theta}z
$$
for some $\theta\in\left[0,2\pi\right)$. Under this transformation, it is straightforward that Eq. $(1)$ becomes
\begin{equation}
e^{-2i\theta}\left(a_1+2ia_2-a_3\right)u_{zz}+2\left(a_1+a_3\right)u_{z\bar{z}}+e^{2i\theta}\left(a_1-2ia_2-a_3\right)u_{\bar{z}\bar{z}}+e^{-i\theta}\left(b_1+ib_2\right)u_z+e^{i\theta}\left(b_1-ib_2\right)u_{\bar{z}}+cu=0.\tag{2}
\end{equation}
Finally, note that the rotational invariance is equivalent to the arbitrariness of $\theta$. Hence compare Eqs. $(1)$ and $(2)$, and the invariance implies the following cases.


*

*If $c\ne 0$, the invariance forces
\begin{align}
a_1+2ia_2-a_3&=0,\\
a_1-2ia_2-a_3&=0,\\
b_1+ib_2&=0,\\
b_1-ib_2&=0.
\end{align}
These results indicate that $a_1=a_3$ and $a_2=b_1=b_2=0$, as is expected.

*If $c=0$ and $a_1+a_3\ne 0$, obviously the same result apply, and we still have the expected conclusion.

*If $c=0$ and $a_1+a_3=0$ and $a_1+2ia_2-a_3\ne 0$ (i.e., $a_1+ia_2\ne 0$), we have
\begin{align}
a_1-2ia_2-a_3&=0,\\
b_1+ib_2&=0,\\
b_1-ib_2&=0,
\end{align}
which, however, do not admit any real solution: the four equalities yield $a_1=a_2=0$ and violate the condition $a_1+ia_2\ne 0$.

*If $c=0$ and $a_1+a_3=0$ and $a_1+2ia_2-a_3=0$, these conditions lead to $a_1=a_2=a_3=0$, making the equation no longer second order.


To sum up, the desired conclusion is completely proven.



*

*Solutions to the Laplace equation and harmonic functions are exactly the same. As you mentioned, one way to define the harmonic functions is to take them as the solutions to the Laplace equation.

*Harmonic functions are not necessarily radial. Radial harmonic functions are called the fundamental solutions to the Laplace equation. In 2-D, it is $\log r$; in 3-D, it is $1/r$. These solutions are essential, and can be made use of to construct Green functions to help solve Poisson equations.

*Let $f$ be a harmonic function, and suppose it yields a separation of variables as $f(r,\theta)=F(r)\Theta(\theta)$. Then $F$ complies with a radial equation, and $\Theta$ is called a spherical harmonic function. These functions are essential in, say, quantum mechanics.

*In general, $f$ can be expressed as
$$
f=\sum_nF_n\Theta_n,
$$
where each $F_n$ complies with a radial equation, and each $\Theta_n$ is a spherical harmonic function. This expression can be obtained by solving the Laplace equation by separation of variables.

