Tensor Field Decomposition in Space time For vector fields in $R^3$ one knows that there exists a unique decomposition of vector fields in to solenoidal (divergence free) and potential parts. A generalization of this Theorem for symmetric tensors of arbitrary order is known for compact Riemannian manifolds with boundary , e.g [Theorem 3.3.2], Sharafutdinov, Integral Geometry of Tensor Fields.
I am curious if such a Theorem is known for any class of Lorentzian Manifolds $M$ too (in particular for Space-time)? i.e. Assume $f$ is a vector field (or more generally a symmetric tensor of order $m$), then does there exist a divergence free vector field (or a symmetric tensor of order $m$) (solenoidal) $f^s$  and a function (or a $m-1$ form) $v$ such that $f=f^s+dv$ , $v|_{\partial M} = 0 $ and where $d$ represents the operator which takes symmetrized co variant derivatives. 
As a first attempt at proof, following Sharafutdinov's approach, we will say that If such a decomposition were to exist then we should be able to find a unique solution to the PDE system $\delta d (v) = \delta f$, $v|_{\partial M} = 0 $. Here $\delta$ represents the divergence operator. Now, if we calculate the symbol for $\delta d$, we can show that at least for the case when $v$ is a function, the operator $\delta d$ is in fact strictly hyperbolic. For higher order tensors, $\delta d$ fails to be strictly hyperbolic. So the question really reduces to the following : When do we know that a solution for a hyperbolic system $Lv=f$, $v|_{\partial M} = 0 $ exists and is unique? 
Also, in the Space-time setting, it looks more natural to work in a non-compact manifold (at least which go in future time directions forever e.g. like a $[0,\infty] \times N$ manifold where $N$ is a compact Riemannian Manifold. Note that Sharafutdinov's theorem is really only for compact Riemannian manifolds, then does it make sense for us to ask a similar question for a non-compact Lorentzian manifolds? Any ideas/comments are welcome as I am feeling out of my depth here!
 A: The class of Lorentzian manifolds that is best suited for this kind of question is that of globally hyperbolic ones. Analyzing the space of solutions of a hyperbolic PDE that does not satisfy an analog of the globally hyperbolic condition is rather messy business.
Let me write in more explicit formulas what you would like to do. For simplicity of notation, let me stick to rank-3 symmetric tensors $f_{abc} = f_{(abc)}$. The generalizations to higher/lower rank tensors should be clear. You are looking to decompose any such tensor as
$$ f_{abc} = 3\nabla_{(a} v_{bc)} + s_{abc} = \nabla_a v_{bc} + \nabla_b v_{ca} + \nabla_c v_{ab} + s_{abc} , $$
where $v_{ab} = v_{(ab)}$, $s_{abc} = s_{(abc)}$, and $\nabla^a s_{abc} = 0$ ($s_{abc}$ is solenoidal). In the Riemannian context, such a decomposition is possible because the composition
$$ 3\nabla^a \nabla_{(a} v_{bc)} = \nabla^a \nabla_a v_{bc} + (3-1)\nabla_{(b} \nabla^a v_{c)a} + (3-1)R^a_{(b} v_{c)a} , $$
with $R_{ab}$ the Ricci tensor, defines an elliptic operator, as is shown in Sharafutdinov. The point is that the second term on the rhs does not destroy the positivity of the symbol of the Laplacian $\Delta = \nabla^a \nabla_a$. As you've noted, Lorentzian situation is more subtle. In fact, it is hard to see how such a second order perturbation of the d'Alembertian $\square = \nabla^a \nabla_a$ defines a hyperbolic operator.
I don't know if you are really constrained to use the above gradient and divergence operators. But in case you are not, I'll propose an alternative one, which will be more convenient. Namely, let us look for a decomposition of the form
$$ f_{abc} = 3\nabla_{(a} v_{bc)} + s_{abc} , $$
where $v_{ab} = v_{(ab)}$, $v_a{}^a = 0$ and $s_{abc} = s_{(abc)}$ such that
$$ (Ds)_{bc} = \nabla^a s_{abc} - \frac{(3-1)}{2}\nabla_{(b} s_{c)a}{}^a = 0.$$
Then, denoting $(\nabla v)_{abc} = 3\nabla_{(a} v_{bc)}$, we have the composition
$$
  (Lv)_{bc} = (D(\nabla v))_{bc} = \square v_{bc} + (3-1) R^a_{(b} v_{c)a} ,
$$
where we have used $(\nabla v)_{ca}{}^a = 2 \nabla^a v_{ca} + (3-2)\nabla_c v_a{}^a = 2 \nabla^a v_{ca}$. Then, any solution $v_{bc}$ of $L[v]_{bc} = (Df)_{bc}$ generate the decomposition $f_{abc} = (\nabla v)_{abc} + s_{abc}$, where $s_{abc} = f_{abc} - (\nabla v)_{abc}$.
Now, the operaor $L$ has the same principal symbol as the d'Alembertian $\square$. Such an operator is called normally hyperbolic. The properties and solvability of normally hyperbolic operators are extensively studied for instance in [BGP] and [W]. Hence, what we can say about the solvability of the equation
$$
  L[v]_{bc} = (Df)_{bc}
$$
is the following. Let us assume that all tensor fields are smooth, but without imposing any support conditions. If $\Sigma$ is a Cauchy surface in the Lorentzian spacetime (the existence of a Cauchy surface is equivalent to the condition of global hyperbolicity) then there is a unique solution $v_{bc}$ to the above equation such that $v_{bc}$ and its first normal derivatives both vanish on $\Sigma$. Further, if $(Df)_{ab}$ has compact or even spacelike compact support, $v_{bc}$ will have spacelike compact support. If you prefer, you can consider instead $\Sigma = \partial \Sigma^+$ as the boundary of the part of the spacetime $\Sigma^+ \cong [0,\infty) \times \Sigma$ to the future of $\Sigma$.
[BGP] Bär, Christian; Ginoux, Nicolas; Pfäffle, Frank MR 2298021 Wave equations on Lorentzian manifolds and quantization, ESI Lectures in Mathematics and Physics ISBN: 978-3-03719-037-1.
[W] Waldmann, Stefan Geometric Wave Equations, [arXiv:1208.4706]
A: Thank you @Igor for your reply yesterday. I re-worked my calculations and here are a couple of questions / comments based on that :
For the case when $v$ is a function : 
$\delta d v = g^{ij}(\partial_j \partial_i v - \Gamma^{k}_{ij} \partial_{k}v)$ 
 which can then be rewritten in the form $(g^{ij} \partial_{i} \partial_{j} v + b^{k}(x) \partial_{k} v)$. This from Michael Taylor, Partial Differential Equations, (8.2) is hyperbolic. In fact it is strictly hyperbolic. In this case, it looks like I should be able to solve the Cauchy problem. 
Assume now that $v$ is a 1-form : 
Now $\delta d v = 1/2 g^{jk}(\partial_k \partial_i v_j + \partial_k \partial_j v_i)$ + lower order things. The principal symbol for this operator is given by $1/2( \langle \xi,\xi \rangle_g Id + \xi(g\xi)^{T})$. On calculating its determinant, one can show that this is not strictly hyperbolic in the sense of Hormander, The Analysis of Linear Partial Differential Operators II, (section 12). In fact, this is the same calculation that @Igor has done. However, from the definition of normal hyperbolic operator as given in [BGP] Wave equations on Lorentzian manifolds and quantization (section 1.5), the operator $\delta d$ in my work seems not to be a normal hyperbolic operator. In such a case, can the PDE system with appropriate Cauchy formulation as suggested by Igor above still be solved?
A: See also 
MR0317200 Barbance, Christiane Sur les tenseurs symétriques. (French) C. R. Acad. Sci. Paris Sér. A-B 276 (1973), A387–A389. 
