Is every gradient vector field a divergence free vector field? What is an example of  a  gradient vector field $X$ on a  Riemannian manifold $(M,g)$ which cannot  be converted to a divergence free vector field via the following processes:

*

*First we remove the singularities $S$ from $M$ then we set $M:=M\setminus S$


*We are allowed  to reparameterize $X$ to $X:=fX$ for some positive function $f$


*We are allowed to change the initial Riemannian metric $g$ to a new metric $g'$ for computation of divergence of $X$ with respect to this new $g'$ to obtain a vector field $X$ on $M$ with $\operatorname{div}_{g'} X=0$.
In other words, with some abuse of terminology, we ask: "Is every gradient vector field a divergence free vector field?"
An obvious example is $X=x\partial_x+y\partial_y$ is a divergence free vector field on the punctured plane after rescaling $X:=\frac{1}{x^2+y^2}X$
Is the answer yes, at least in low dimensions?
 A: I think the answer is no as soon as your gradient vector field admits a saddle point where the divergence is non-zero.
Let $\omega$ denote the volume form associated to the Riemann metric. We have
$$\mathrm{div}(X) \omega = X\cdot \omega$$
where $X\cdot \omega$ denotes the Lie derivative. The goal is to find positive functions $f$ and $g$ such that
$$(fX)\cdot (g\omega) = X\cdot(fg)~\omega + (fg) \mathrm{div}(X)~\omega = 0~.$$
In other words, we want the function $h=\log(fg)$ to satisfy
$$X\cdot h = -\mathrm{div}(X)~.$$
This is a dynamical question: we ask whether the function $\mathrm{div}(X)$ is a coboundary along the flow of $X$. Of course a gradient flow does not have very rich dynamics, but a saddle point is already too much for the following reason:
Assume $X$ has a saddle point. Then one can find sequences $(x_i)$ and $(y_i)$ which are bounded in $M\backslash S$ such that $y_i$ is on the trajectory of $x_i$  along the flow of $X$, and such that the trajectory from $x_i$ to $y_i$ is very long and spends most of its time very close to the saddle point $s$.
Assume now that we have $h:M\backslash S \to \mathbb R$ such that $\mathrm{div}(X)= X\cdot h$. Then
$$\int_{x_i}^{y_i} \mathrm{div}(X) = h(y_i)-h(x_i)$$
is bounded independently of $i$. (Here $\int_{x_i}^{y_i} \mathrm{div}(X)$ denotes the integral of the function $\mathrm{div}(X)$ along the trajectory of $X$ from $x_i$ to $y_i$.)
But, on the other side, since this trajectory spends a very long time close to $s$, we have that
$$\int_{x_i}^{y_i} \mathrm{div}(X)\underset{i\to +\infty}{\longrightarrow} \pm \infty$$
as soon as $\mathrm{div}(X)(s) \neq 0$. This is a contradiction.
