Existence of second order potential for PDE 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)?
 A: First of all, you have a sign wrong in your formula for the curvature.  The curvature tensor you gave has positive constant sectional curvature +1 while you claim that you want negative sectional curvature (i.e., hyperbolic space), which would flip the sign of $R$.  Second, when you speak of spherical harmonics, I believe you must be copying the formula for the unit $n$-sphere, $S^n$, not hyperbolic space $H^n$.  This also causes an error in your potential formula (2), which would be correct for the $n$-sphere, but should be $T_{ab} = D_aD_bT - T h_{ab}$ for hyperbolic space.
I'll give the answer for hyperbolic space, since that is what you want, but be aware that you'll need to flip signs to get the same answer for the $n$-sphere.
The result you are seeking follows immediately from the Frobenius theorem, using the techniques mentioned in the MO question you cite at the end.  The idea is simply this:  Take the tensor $T$ as given.  Let $\omega_i$ be any $h$-orthonormal frame field (which can be chosen globally on $H^n$ since it is contractible) and let $\theta_{ij}=-\theta_{ji}$ be the unique $1$-forms satisfying $\mathrm{d}\omega_i = -\theta_{ij}\wedge\omega_j$. (Here and below, I am using the 'Einstein' summation convention.) By the assumption that the sectional curvature is identically $-1$, we have $\mathrm{d}\theta_{ij} = -\theta_{ik}\wedge\theta_{kj} - \omega_i \wedge\omega_j$.
Now, on $X = H^n\times\mathbb{R}\times\mathbb{R}^n$, with projections $u:X\to\mathbb{R}$ and $(u_i):X\to\mathbb{R}^n$ onto the second and third factors, consider the Pfaffian system $\mathcal{I}$ generated by the $(n{+}1)$ linearly independent $1$-forms
$$
\xi = \mathrm{d}u - u_i\ \omega_i
\quad\text{and}\quad
\xi_i = \mathrm{d}u_i +\theta_{ij}\ u_j - (T_{ij}+u\,\delta_{ij})\ \omega_j\,.
$$
By the hypotheses on $T = T_{ij}\omega_i\omega_j$ and the curvature of $h$, this system is Frobenius, i.e., it satisfies $\mathrm{d}\xi \equiv \mathrm{d}\xi_i\equiv 0 \mod \mathcal{I}$.  
Thus, $X$ is foliated by the leaves of $\mathcal{I}$, which are transverse to the fibers of the projection $\pi:X\to H^n$ onto the first factor.  Since the system is affine linear in $(u,u_i)$, it follows that each leaf $L\subset X$ of $\mathcal{I}$ becomes a covering space of $H^n$ under the projection $\pi:L\to H^n$.  Since $H^n$ is connected and simply-connected, such a projection is a diffeomorphism of the leaf $L$ with $H^n$ and hence has an inverse, which can be written in the form $\sigma:H^n\to L\subset X$ of the form $\sigma(p) = \bigl(p,f(p),f_i(p)\bigr)$.  By construction, the function $f = u\circ\sigma:H^n\to\mathbb{R}$ and $(f_i) = (u_i)\circ\sigma:H^n\to\mathbb{R}^n$ must satisfy
$$
df = f_i\ \omega_i
\qquad\text{and}\qquad
df_i = -\theta_{ij}\ f_j + (T_{ij}+f\delta_{ij})\ \omega_j\,.
$$
The function $f$ is the potential that you seek.
