Proof that $a\nabla^2 u = bu$ is the only homogenous second order 2D PDE unchanged/invariant by rotation

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.

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.

---

1. 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.
2. 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.
3. 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.
4. 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.