Given some function $f \in L^2(\mathbb{R}^2)$, I'm interested in finding a positive semi-definite differential operator $\mathcal P: L^2(\mathbb{R}^2) \rightarrow L^2(\mathbb{R}^2)$ that is quadratic in $f$ and invariant under the the action of $\textrm{SL}{(2, \mathbb{R})},$ such that $\forall A \in \textrm{SL}(2, \mathbb{R})$ and $\forall {\bf x} \in \mathbb{R}^2,$ $$ {\mathcal P} f(A {\bf x}) = [{\mathcal P} f] (A {\bf x} ).$$
After thinking for some time, I've come up with two operators that are invariant and P.S.D, but not quadratic.
For example, suppose we consider the operator $${\mathcal P} = \left( \frac{\partial^2}{\partial x^2} \frac{\partial^2}{\partial y^2} - \left[\frac{\partial^2}{\partial x y}\right] \right)^2,$$ which is the squared determinant of the Hessian. It's clear that it is P.S.D. and invariant under transformations in $\textrm{SL}(2, \mathbb{R})$, though it is quartic in $f$.
Furthermore, letting $H$ denote the Hessian and $J \in \textrm{SO}(2)$ be a rotation by $90^\circ$, the operator $$ {\mathcal P} = \left(\nabla^T J^T H \ J \ \nabla\right)^2,$$ is also invariant and P.S.D., but is not quadratic.
I'm asking this question in the hope that someone might know of such a quadratic P.S.D. differential operator that is invariant under $\textrm{SL}(2, \mathbb{R})$ (if it even exists) or be able to point me toward a with a few other ideas I could try.
Some possibly related question(s):
Projective invariant differential operator
Classification of diffeomorphisms by association of differentials with Lie groups
I will take your question seriously but not literally since you are asking for a differential operator $D: L^2({\mathbb R}^2)\to L^2({\mathbb R}^2)$ and there are no differential operators $D$ of order $>0$ which take all $L^2$-functions to $L^2$-functions: You would need distributions as values of $D$.
Thus, I will assume that $L^2({\mathbb R}^2)$ in your question means $C^\infty({\mathbb R}^2)$. Then you get your example:
Take $Z=x\frac{\partial}{\partial x} + y \frac{\partial}{\partial y}$: As a vector-field, it sends each point with coordinates $(x,y)$ to vector with coordinates $(x,y)$, this is why $Z$ is invariant under the action of $GL(2, {\mathbb R})$. Then take $D=Z\otimes Z$. As a differential operator, it acts on smooth functions by $$ D: f\mapsto (x\frac{\partial f}{\partial x} + y \frac{\partial f}{\partial y})^2. $$
One can prove that among strictly 1st order strictly quadratic PSD differential operators, up to scalar, this is the only one which is $SL(2, {\mathbb R})$-invariant.
The same works in higher dimensions as well, your $GL(n, {\mathbb R})$-invariant differential operator will be $$ Z\otimes Z, Z=\sum_{i=1}^n x_i\frac{\partial}{\partial x_i}. $$
If I were to take your question literally but not seriously, my answer would be $$ D: L^2({\mathbb R}^2)\to L^2({\mathbb R}^2), D(f)=a f^2 $$ where $a\ge 0$ is a fixed constant. Such $D$ is a PSD, continuous, quadratic differential operator of order 0.