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 and consider the system of differential equations relative to any $h$-orthonormal frame field $\omega_i$ (which can be defined globally since the hyperboloid is diffeomorphic to $\mathbb{R}^n$): $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + (T_{ij}+f\delta_{ij})\ \omega_j\,, \tag1 $$ where $\theta_{ij}=-\theta_{ji}$ are the unique $1$-forms satisfying $\mathrm{d}\omega_i = -\theta_{ij}\wedge\omega_j$. This is an affine system of total differential equations for the function $f$ and its derivatives $f_i$ relative to the frame field. The hypotheses on $T_{ij}$ imply that this system is Frobenius and hence has local solutions everywhere. The homogeneous system (corresponding to $T=0$) is $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + f\delta_{ij}\ \omega_j\,, \tag2 $$ and this has an $n{+}1$-dimensional space of solutions on any connected open set. (Basically, these are the functions used to embed hyperbolic space into Lorentzian $(n{+}1)$-space as a space-like hyperquadric.)

Because $H^n$ is simply connected and the homogeneous system is linear, the usual cohomological argument shows that the system $(1)$ always has a global solution $f$, which is unique up to the addition of a solution of $(2)$. Such a global solution $f$ is the potential $T$ that you seek.

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.

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, 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 and consider the system of differential equations relative to any $h$-orthonormal frame field $\omega_i$ (which can be defined globally since the hyperboloid is diffeomorphic to $\mathbb{R}^n$): $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + (T_{ij}+f\delta_{ij})\ \omega_j\,, \tag1 $$ where $\theta_{ij}=-\theta_{ji}$ are the unique $1$-forms satisfying $\mathrm{d}\omega_i = -\theta_{ij}\wedge\omega_j$. This is an affine system of total differential equations for the function $f$ and its derivatives $f_i$ relative to the frame field. The hypotheses on $T_{ij}$ imply that this system is Frobenius and hence has local solutions everywhere. The homogeneous system (corresponding to $T=0$) is $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + f\delta_{ij}\ \omega_j\,, \tag2 $$ and this has an $n{+}1$-dimensional space of solutions on any connected open set. (Basically, these are the functions used to embed hyperbolic space into Lorentzian $(n{+}1)$-space as a space-like hyperquadric.)

Because $H^n$ is simply connected and the homogeneous system is linear, the usual cohomological argument shows that the system $(1)$ always has a global solution $f$, which is unique up to the addition of a solution of $(2)$. Such a global solution $f$ is the potential $T$ that you seek.

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 and consider the system of differential equations relative to any $h$-orthonormal frame field $\omega_i$ (which can be defined globally since the hyperboloid is diffeomorphic to $\mathbb{R}^n$): $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + (T_{ij}+f\delta_{ij})\ \omega_j\,, \tag1 $$ where $\theta_{ij}=-\theta_{ji}$ are the unique $1$-forms satisfying $\mathrm{d}\omega_i = -\theta_{ij}\wedge\omega_j$. This is an affine system of total differential equations for the function $f$ and its derivatives $f_i$ relative to the frame field. The hypotheses on $T_{ij}$ imply that this system is Frobenius and hence has local solutions everywhere. The homogeneous system (corresponding to $T=0$) is $$ df = f_i\ \omega_i \qquad\text{and}\qquad df_i = -\theta_{ij}\ f_j + f\delta_{ij}\ \omega_j\,, \tag2 $$ and this has an $n{+}1$-dimensional space of solutions on any connected open set. (Basically, these are the functions used to embed hyperbolic space into Lorentzian $(n{+}1)$-space as a space-like hyperquadric.)

Because $H^n$ is simply connected and the homogeneous system is linear, the usual cohomological argument shows that the system $(1)$ always has a global solution $f$, which is unique up to the addition of a solution of $(2)$. Such a global solution $f$ is the potential $T$ that you seek.