Conformal transformation of the divergence Let $f$ be a smooth function on a $n$-dimensional Riemannian mainfold $(M, g)$, so that $\tilde{g} = e^{2f} g$ is a conformal transformation of $g$. I'm trying to show that the divergence operator transforms as
$$
\text{div}_{\tilde{g}} X = \text{div}_{g} X + n e^{-2f} X(f)
$$
but I'm a little stuck. Dropping down to coordinates, the divergence looks like
$$
\text{div}_\tilde{g} X = \sum_i \frac{1}{\sqrt{\det \tilde{g}}} \frac{\partial}{\partial x^i} \left( X^i \sqrt{\det \tilde{g}} \right) = \sum_i \frac{e^{-nf}}{\sqrt{\det g}} \frac{\partial}{\partial x^i} \left( X^i e^{nf} \sqrt{\det g} \right)
$$
An application of the product rule to the differential term gets you one term looking like $\text{div}_g X$ (hooray!) and one term looking like $nX(f)$. This isn't what I want - it's off by the factor of $e^{-2f}$. What's going wrong with this very simple calculation? I'm pretty sure the formula at the top is correct.

Also, is there a coordinate-free way to do this? Maybe that would help.
Here's an relatively straightfoward, coordinate-free way to approach the problem:

The divergence of a metric $g$ actually depends only on its volume form $\mu$ (for either orientation, it doesn't matter which). It can be characterized by
$$\mathcal{L}_X \mu = (\text{div } X) \mu.$$

Now, the volume form $\widehat{\mu}$ of a conformally related metric $\widehat{g} = e^{2f} g$ is
$$\widehat{\mu} = e^{n f} \mu;$$
substituting in the above characterization of the divergence---and using the Leibniz identity $$\mathcal L_X (h \alpha) = (X \cdot h) \alpha + h \mathcal L_X \alpha$$ for a function $h$ and form $\alpha$---gives an identity relating the divergences of the two metrics.
Your computation is correct but I think what you want to show is incorrect: divergence does not change if you do a simple rescale $g \to c^2g$ when $c$ is constant, but the term $ne^{-2f}X(f)$ does change if you replace $f$ by $f+C$ (that is one more simple rescale after the conformal change).