For fixed $\beta$ in the sphere $\mathbb S^2\subset \mathbb R^3$, we can define a "zonal derivative" via the distribution $\delta'(\langle \cdot, \beta \rangle)$ on $\mathbb S^2$ -- here $\delta'$ is the derivative of the Delta function, and $\langle \cdot, \cdot \rangle$ is the usual Euclidean scalar product on $\mathbb R^3$.
I'd like to combine two such distributions. To be more precise, given test function $u \in C^\infty(\mathbb S^2)$, I am interested in knowing more about the distribution that maps $u$ to
$ \tilde u(x) = \int_{\mathbb S^2} \int_{\mathbb S^2}\hspace{0.5em} u(\beta)\,\,\, \delta'(\langle \alpha, \beta \rangle) \,\,\, \delta'(\langle x, \beta \rangle) \,\,\,\,dS(\alpha) dS(\beta) .$
Here $dS$ is the (Euclidean) surface measure on the sphere.
My question: what can be said about $u \mapsto \tilde u$ as a differential operator? What order does this distribution have? This kind of object shows up very naturally in certain x-ray reconstruction problems, in particular in the more general form with weighting function $a \in C^\infty(\mathbb S^2 \times \mathbb S^2 )$:
$ \tilde u_a(x) = \int_{\mathbb S^2} \int_{\mathbb S^2}\,\,\,a(\alpha,x)\,\, u(\beta) \,\,\delta'(\langle \alpha, \beta \rangle) \delta'(\langle x, \beta \rangle)\,\,\, dS(\alpha) dS(\beta) .$
The tomographic reconstruction papers I've read using something like this usually go into coordinates quite quickly and are written from a very applied point of view. I'd be very interested in finding out more from the distribution-theoretic side.