First order PDE system on a complete Riemannian manifold Let $X_1, X_2$ be two orthogonal everywhere non-zero vector fields on a complete Riemannian manifold $M$. Can one always solve the solve the system $$X \phi = 1, Y\phi = 0,$$ 
$$X\psi = 0, Y\psi = 1 ?$$
where $\phi, \psi \in C_0^\infty (M)$? I am not too sure this can be done always, but may be there are good sufficient conditions known for such linear systems?
 A: I assume that you're using $C^\infty_0(M)$ with its usual meaning of the set of smooth, compactly supported, real-valued functions on $M$. If so, then the answer is that there are never global solutions.
If $M$ is noncompact and $\phi\in C^\infty_0(M)$, then there will be an open subset of $M$ where $\phi\equiv 0$, which contradicts $X\phi = 1$. On the other hand, if $M$ is compact, then $X$ and $Y$ are complete, meaning that every integral curve is defined for all $t\in \mathbb R$. If $\gamma\colon \mathbb R\to M$ is an integral curve of $X$, say, then 
$$
\frac{d}{dt} \phi(\gamma(t)) = X\phi(\gamma(t)) = 1,
$$
so $\phi(\gamma(t)) = t+c$ for some constant $c$. This contradicts the fact that every continuous function on a compact space is bounded.
So there are no compactly supported solutions. Suppose we look instead for global solutions that are not necessarily compactly supported. If $M$ is compact, every function is compactly supported, so the argument above shows that there are no global solutions in that case. If $M$ is noncompact, there still might be no global solutions. For example, as @Anthony Carapetis mentioned, there will be no solutions if $X$ or $Y$ has any closed integral curves. In fact, if $X$ has an integral curve that remains in a compact subset of $M$ for all positive or negative time, then the same argument leads to a contradiction, so again there are no solutions. The question of existence of global solutions is thus subtly tied up with the topology of $M$ and of the foliations determined by $X$ and $Y$, and I doubt that there's any simple general answer.
So that brings us to the question of local solutions. As @Anthony Carapetis mentioned, a sufficient condition for the existence of local solutions in a neighborhood of each point is $[X,Y]\equiv 0$. One reference for this is Theorem 9.46 in my Introduction to Smooth Manifolds. 
When the dimension of $M$ is $2$, the same theorem shows that this is also necessary, because any such pair of functions $(\phi,\psi)$ forms local coordinates on $M$. (And, incidentally, it should be noted that this implies the given Riemannian metric is flat, so flatness of the metric is another necessary condition. See Theorem 13.14 in the same book.)
However, in higher dimensions, the condition $[X,Y]\equiv 0$ is not necessary. Here's an example in $3$ dimensions for which there exist global solutions. Define vector fields $X$ and $Y$ on $\mathbb R^3$ by
$$
X = \frac{\partial}{\partial x} - y\frac{\partial}{\partial z},
\qquad 
Y = \frac{\partial}{\partial y} + x\frac{\partial}{\partial z}.
$$
Then $[X,Y] \ne 0$, but $\phi(x,y,z) = x$ and $\psi(x,y,z)=y$ are global solutions to the problem. (If you want to relate this to a Riemannian metric, just let $g$ be the metric for which $\{X,Y,\partial/\partial z\}$ is an orthonormal basis.)
So one might wonder if (local) solutions always exist in higher dimensions. But that's not the case either: For example, again in $\mathbb R^3$, if we set 
$$ 
X = \frac{\partial}{\partial x},
\qquad 
Y = \frac{\partial}{\partial y} + x\frac{\partial}{\partial x},
$$
then $[X,Y]=X$, and any solution $\phi$ would have to satisfy 
$$
1 = X\phi = [X,Y]\phi = X(Y\phi) - Y(X\phi) = X(0) - Y(1) = 0,
$$ 
a contradiction.
I haven't figured out whether there's a simple necessary and sufficient condition for the existence of local solutions in higher dimensions. Interesting question.
