Quadratic P.S.D. differential operator that is invariant under $\textrm{SL}(2, \mathbb{R})$ 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
Proof that $a\nabla u = b u$ is the only homogenous second order 2D PDE unchanged/invariant by rotation
 A: 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.
