In coordinate-free language, my question is as follows. Let $M$ be an $n$-dimensional manifold with volume form, and let $\mathcal D$ be a smooth (integrable, if necessary) distribution with constant rank $k\leq m$ (by which I mean that $\mathcal D$ is a subbundle of the tangent bundle, spanned by $k$ many every-linearly-independent smooth vector fields). Given a point $p\in M$, is there necessarily a neighborhood $U\ni p$ and vector fields $w_1,\dots,w_k$ on $U$ so that each $w_a$ is divergence-free (with respect to the volume form) and so that $\mathcal D|_U = \text{span}\{w_1,\dots,w_k\}$?

Since my question is local, I will ask it again on $\mathbb R^n$ in coordinates. I will let $x^i$, $i=1,\dots,n$ be the standard coordinates on $\mathbb R^n$. Suppose you are given smooth vector fields $v_a(x) = \sum_{i=1}^n v_a^i(x) \frac{\partial}{\partial x^i}$ for $a=1,\dots,k$, where $k\leq n$, and suppose moreover that for each $x$, the set of vectors $\{v_1(x),\dots,v_k(x)\}$ is linear independent. (Put another way, $v$ is an everywhere-full-rank $(k\times n)$-matrix-valued function.) I'm looking for a smooth ${\rm GL}(k)$-valued-function $m(x) = \{m^a_b(x)\}_{a,b=1}^k$ in a neighborhood of $0\in \mathbb R^n$ so that for each $b=1,\dots,k$, $$ \sum_{i=1}^n \sum_{a=1}^k \frac{\partial}{\partial x^i}\bigl[ m^a_b(x)\,v^i_a(x) \bigr] = 0 $$

If such an $m$ (or, if you prefer, $\{w_1,\dots,w_k\}$) exists, the natural questions are:

- Locally, how unique is it? Certainly it can by multiplied by any constant invertible $k\times k$ matrix. Is this it? I.e., if we fix the matrix (spanning set) at $0$ ($p$), does it fix the value in an open neighborhood?
- Can "local" be replaced by "global"? (I'm not even sure whether a rank-$k$ smooth distribution, which for me is a subbundle that's locally spanned by $k$ everywhere-independent vector fields, can be globally spanned by $k$ everywhere-independent vector fields.)

Finally, I should remark that for my application, my distribution is integrable — does this matter for whether the answer is yes? Also, for my application, not only do I only need the result locally, but I actually am working in the formal neighborhood of a point, so for my purposes it's fine if all smooth functions are replaced by formal power series. Occasionally, but rarely, questions of existence of solutions to differential equations are easy in formal land but hard even in smooth land.

If you're feeling extra smart, the last generalization I'd be interested in knowing the answer to is what happens if the rank of the distribution can drop. I.e. fix $k\leq n$, and suppose that we have vector fields $v_1,\dots,v_k$, which may be linearly dependent at $0$. When is there an (invertible) $(k\times k)$-matrix-valued function $m$ so that the above holds? The answer is a resounding "no" when $k=n=1$, but I have no intuition for higher values.

My motivation is that integrable distributions are the most general thing (that I know of, anyway), that can act by smooth "symmetries" of a manifold. I would like to know if every smooth symmetry can be taken to be volume-preserving. For a different interpretation of this general question, see e.g. Moser, On the volume elements on a manifold, Trans. Amer. Math. Soc., 1965, vol. 120, pp. 286--294, MR0182927, in which it is proven that any globally-volume-preserving diffeomorphism is smoothly isotopic through globally-volume-preserving diffeomorphisms to a locally-volume-preserving diffeomorphism.