Does a (smooth, constant-rank, integrable) distribution have a (local) basis of divergence-free vector fields? 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.
 A: For any non-zero vector field $v$ and a volume form $\omega$ there exists (locally) a positive function $f$ such that $L_{fv}\omega=0$ (Indeed, one can take coordinates such that $v=\frac{\partial}{\partial x_1}$, $\omega=A(x_1,...,x_n)dx_1\wedge ...\wedge dx_n$. In these coordinates the condition $L_{fv}\omega=0$ is equivalent to an ordinary differential equation $\frac{\partial f}{\partial x_1}+\frac{1}{A}\frac{\partial A}{\partial x_1}f=0$, I am using the homotopy formula $L_u=i_ud+di_u$ here with $u=fv$). Hence, for any distribution (of constant rank) there exists a local basis of divergent-free vector fields - take any basis of vector fields and multiply each vector field from it on a suitable $f$.
A: As already pointed out in the comments, one does not necessarily have nowhere vanishing vector fields tangent to a foliation.
Although in the differentiable setting this is always the case for orientable one dimensional foliations. But for a foliation $\mathcal F$ to be defined by a divergence free vector field is a very strong condition. In this case the contraction of the volume form with the vector field is a closed $(n-1)$-form, which we will call $\omega$. If we take a leaf $L$ of $\mathcal F$,  two germs of transversals to it, and a path on $L$
connecting the two transversals then the holonomy  map will preserve the measure induced on the transversals by the closed $(n-1)$-form $\omega$. This is a simple application of Stokes Theorem.
In particular, if you have a vector field on $\mathbb R^2$ having a periodic orbit accumulated by non closed orbits then it is impossible to define the underlying  foliation by a divergence vector field on any neighborhood of the periodic orbit. containing the limit cycle
