Does for every vector field there always exist a volume form for which the vector field is a homothety?

If there is such a volume form $\Omega$, by Moser's theorem we can pick local coordinates in which $\Omega=dx^1 \wedge \dots \wedge dx^n$. If in some coordinates $v$ vanishes to order $k$, for some $k>1$, but not at order $k$, then the same is true in any coordinates. For example, we can suppose that $v=f^2w$ with $f=0$ at some point where $df \ne 0$ and $w \ne 0$ and $w$ not tangent to $f=0$ at that point. This condition is coordinate invariant, and we calculate that $\mathcal{L}_v \Omega=C \Omega$ just when $C=0$ and $2\mathcal{L}_w f + f \, \text{div} w = 0$. But at $f=0$ this forces $\mathcal{L}_w f =0$ tangent, a contradiction. So $v$ does not preserve $\Omega$. But then $v$ does not preserve any volume form.

There may be obstructions, even in the projective case and even when the vector field does not vanish to second order. Suppose given a projective structure on an $n$-manifold $M$ and a projective vector field $X$ that vanishes at a point $p\in M$. Then, in projective normal coordinates (either the Thomas-Veblen version or the Cartan version) $x = (x^i)$ centered on $p$, the vector field $X$ will have the form $$ X = a^i_j\ x^j\ \frac{\partial\ \ }{\partial x^i} + (b_j x^j)\ \left(x^i\ \frac{\partial\ \ }{\partial x^i}\right) $$ for some constants $a^i_j$ and $b_j$. If $\Omega = f\,dx^1\wedge\cdots\wedge dx^n$ is an $n$-form on a neighborhood of $p$, then the condition that ${\mathcal{L}}_X\Omega = C\ \Omega$ is that $f$ satisfy the partial differential equation $$ a^i_j\ x^j\ \frac{\partial f}{\partial x^i} + (b_j x^j)\ \left(x^i\ \frac{\partial f}{\partial x^i}\right) +\bigl( (n{+}1) (b_j x^j) + a^i_i - C\bigr) f = 0. $$ Evaluating this at $x=0$, gives $(a^i_i-C)f(0) = 0$, so if $\Omega$ is to be nonvanishing at $0$, one must have $C = a^i_i$. Supposing this, the equation simplifies to $$ a^i_j\ x^j\ \frac{\partial f}{\partial x^i} + (b_j x^j)\ \left(x^i\ \frac{\partial f}{\partial x^i}\right) +(n{+}1) (b_j x^j) f = 0. $$ Writing $f = f_0 + f_1 + \cdots$, where $f_k$ is the $k$-th homogeneous term in the Taylor series, the above equation now implies the recursive relations $$ a^i_j\ x^j\ \frac{\partial f_k}{\partial x^i} + (n{+}k) (b_j x^j) f_{k-1} = 0 \qquad (k\ge 1). $$ Obviously, this cannot be satisfied for $k=1$ with $f_0\not=0$ unless $b_j = a^i_jc_i$ for some constants $c_i$. (Note, in particular, that a vector field $X$ of the above form that does not satisfy this condition provides a negative answer to Vladimir's remaining question; the divergence at $p$ is not really relevant.) Conversely, if this condition holds, then $$ f = (1 + c_ix^i)^{-(n+1)} $$ satisfies the equation, so the desired $n$-form $\Omega$ does exist, with $C = a^i_i$, at least in an open neighborhood of $p$.

Thus, following on Vladimir's comments above, the answer to the motivating question from projective geometry comes down to whether, when the projective curvature does not vanish at $p$, one always has $b_j=a^i_jc_i$ for some constants $c_i$ when $X$ is a projective vector field vanishing at $p$. (Perhaps Vladimir knows whether this is true. Vladimir?)

For a vectorfield $v$ in a region of $\mathbb R^n$ and a volume element $\Omega=f\,\Lambda$ (where $\Lambda$ is the Lebesgue volume element and $f$ is a strictly positive smooth function), the equation $$\mathcal L_v\Omega=C\,\Omega$$ (with $C\in\mathbb R$ constant) is equivalent to the fact that for each measurable set $A$ we have $$|v^t(A)|_\Omega=\mathrm e^{C\,t}|A|_\Omega,$$ where $v^t$ is the diffeomorphism obtained by following the vectorfield $v$ during a time $t$, and $|A|_\Omega:=\int_A\Omega$ is the measure of $A$ according to $\Omega$. Observe that in this formulation we may allow more general density functions $f$.

Consider near $0\in\mathbb R^2$ the vectorfield $$v(x,y)=y\,(x\,\partial_y-y\,\partial x).$$ This vectorfield flows along the semicircles $$\sqrt{x^2+y^2}=const,\quad y­­>0$$ towards the left, from a geometrically repelling fixedpoint on the positive $x$-axis to a similar attracting fixedpoint on the negative $x$-axis. Each semiannulus $$\epsilon\leq\sqrt{x^2+y^2}\leq 2\epsilon,\quad y\geq 0$$ is preserved, but its content is pressed towards the left (at an exponential rate). We see immediately that the volume element $\Omega$ cannot exist for this $v$, because the right half of the semiannulus is being expanded and the left half is being contracted. But this vectorfield still has zero divergence at the origin.

To get an example with nonzero divergence we go to $\mathbb R^3$ and define $$v(x,y,z)=y\,(x\,\partial_y-y\,\partial x)-z\,\partial_z,$$ which has nonzero divergence in a neighbourhood of the origin. The flow of this vectorfield sends each semi-solidtorus $$S_k=\left\{(x,y,z):\epsilon\leq\sqrt{x^2+y^2}\leq 2\epsilon,\ y\geq 0,\ z\in\left[\frac \epsilon{\mathrm e^{k+1}},\frac \epsilon{\mathrm e^k}\right]\right\}$$ to the following semi-solidtorus $S_{k+1}$ after one unit of time. But the content is again compressed towards the left. So the density function $f$ has singularities on the plane $z=0$.

In more detail, by comparing the Lebesgue measures $|S_k|_\Lambda=c\,e^{-k}$ with the $\Omega$-measure $|S_k|_\Omega=c'\,C^k$ we see that we need to have $C=\frac 1{\mathrm e}$ to hold a hope that the density factor $f$ is bounded above and below by strictly positive constants. But this value is actually not important because we needn't care know how much each $S_k$ measures according to $\Omega$. It's enough to know that the density on the left is much (exponentially in $k$) greater than the density on the right. We will have (exponential) singularities of $f$ along the positive $x$-axis or the negative $x$-axis.