Let me start with a definition: Given a Riemannian manifold (M,g) [oriented, connected] we call a smooth vector field $Y$ concircular with potential function $f\in C^{\infty}(M)$ if for every smooth vector field $X$ we have the relation $\nabla_XY=f X$ with $\nabla$ denoting the Levi-Civita connection.
Examples are the position vector field on $\mathbb{R}^n$ with $\nabla_X\vec{x}=X$ and the standard basis $e_i$, $\nabla_Xe_i=0$.
Note that the $e_i$ are Killing (meaning their local flow preserves the metric) and that the corresponding potential function is identically zero. It is not hard to show that this is true in general: If $Y$ is a smooth Killing field on some Riemannian manifold, which is also concircular, then its potential function must be identically zero.
Now Killing fields in particular preserve the volume. So I was wondering what happens if we relax the assumption and consider divergence-free vector fields $Y$ which are concircular. Does the potential function then vanish necessarily as well?
Or are there examples of divergence-free concircular vector fields with non-vanishing potential function?
I tried to plug in $X=Y$ in the defining equation and apply the divergence on both sides and use some calculus formulae, but that didn't lead me anywhere.
I am happy with results for the 3-d case, in case the dimension is relevant.
By now I found an answer. So this is for everyone who's interested.
In the following $(M,g)$ is an arbitrary oriented, connected smooth Riemannian $n$-manifold. We assume that $Y$ is a smooth, concircular vector field.
It is easiest to work in normal coordinates. Fix any $p\in M$ and a normal coordinate system centered around $p$. We denote the induced tangent basis vectors by $\partial_i$. So fix any $1\leq i \leq n$ and set $X=\partial_i$ (strictly speaking this is only defined on some open neighbourhood, but we can multiply it by an appropriate bump function around $p$ to get a vector field defined on all of $M$). Then due to the conciruclar condition we find:
$(\nabla_{\partial_i}Y)(p)=f(p)\partial_i(p)$.
Now multiply the above identity by $\partial_k(p)$ for some fixed $1\leq k\leq n$
$g(\partial_k(p),(\nabla_{\partial_i}Y)(p))=f(p)\delta_{ik}$,
where we used that the $\partial_i$ are orthonormal at $p$. Now locally in normal coordinates the left hand side reads: $(\partial_iY^k)(p)$. Thus for $k\neq i$, we have $(\partial_iY^k)(p)=0$. This in particular implies that the associated $1$-form $\alpha_{Y}$ (via the Riemannian metric) is closed. In particular in $3$-d, every concircular vector field must be irrotational. Now if we set $k=i$ and take the sum over all $i$, we obtain (using the summation convention)
$(\partial_iY^i)(p)=nf(p)$
The left hand side is now precisely the divergence of $Y$ expressed in normal coordinates. Hence we arrive at
$\text{div}(Y)=nf$.
So a concircular vector field is divergence-free if and only if its potential function is identically zero.
This then in turn implies the following equivlances
Let $(M,g)$ be as above and $Y$ a smooth, concircular vector field.
$Y$ is divergence-free $\Leftrightarrow$ The potential function is indetically zero $\Leftrightarrow$ $Y$ is a Killing vector field $\Leftrightarrow$ $Y$ is a parallel vector field.