There is a statement in the literature (see the paragraph between equations (18) and (19) in http://aip.scitation.org/doi/10.1063/1.523863), which I would like to generalise, but I don't have a nice proof of the original claim. The statement is the following:

Given any symmetric tensor field $T_{ab}$ on a 3-dimensional hyperboloid $H$ with \begin{equation} D_{[a}T_{b]c}=0, \tag{1} \end{equation} there exists a scalar field $T$ on $H$ satisfying \begin{equation} T_{ab} = D_a D_b T + T h_{ab}, \tag{2} \end{equation} where $h_{ab}$ is the metric on the hyperboloid and $D$ is the covariant derivative on the hyperboloid. The Riemann tensor on $H$ satisfies $R_{abcd} = h_{ac}h_{bd} - h_{ad} h_{bc}$.

(I am also happy to assume that $T_{ab}$ is traceless, which together with (1) implies that $T_{ab}$ is divergence-less.)

Physically this is saying the the tensor $T_{ab}$ admits a (second order) potential $T$. I would like to generalise it to other dimensions, possibly to other manifolds, etc. But the only proof of this statement which I know, relies on explicit use of spherical harmonics on $H$.

Are there any general methods to determine if given tensor satisfying certain PDE can be written as some differential operator (possibly involving geometric data) acting on some other tensor?

After some search I found Characterizing Hessians among symmetric bilinear tensors, which I can mimic to show that given $T$ then $T_{ab}$ constructed in (2) indeed solves (1). But how can I show that all solutions to (1) are given by (2)?