Formula for divergence of conformal metric If we have riemannian metrics $g, \tilde{g}$ on a 3-manifold $M$ and a differentiable function $\psi$ such that $g_{ij} = \psi^4 \tilde{g}_{ij}$, then for any symmetric and traceless tensor field $A^{ij}$ it holds that
$$D_jA^{ij} = \psi^{-10}\tilde{D}_j(\psi^{10}A^{ij}).$$
$D$ denotes the covariant derivatives. This identity should easily follow from known identities for conformal metrics, but my computations just don't add up. Does anyone know how to do the calculations?
 A: Let us change the notation slightly, so that $\widehat{g} = \Omega^2 g$ is the rescaled Riemannian metric, where $\Omega = e^{\omega}$ and $\omega \in C^{\infty}(M)$. In particular, $\Omega > 0$ everywhere on $M$. For the meantime, let us work in a greater generality, assuming that $M$ is a closed smooth manifold of dimension $n = \dim M$. As in the question, $g$ is a Riemannian metric on $M$.
Extending this notation, we write $\nabla$ and $\widehat{\nabla}$ for the Levi-Civita connections of the metrics $g$ and $\widehat{g}$, respectively.
As we know,
$$
\widehat{\nabla}_X Y = \nabla_X Y + (X \omega )Y + (Y \omega )X - g(X,Y) \operatorname{grad}\omega \tag{1}
$$
Using the abstract index notation, we can rewrite the above equation in the following form:
$$
X^a \widehat{\nabla}_a Y^b = X^a \nabla_a Y^b + (X^a \nabla_a \omega) Y^b + (Y^a \nabla_a \omega) - X^a Y^c g_{a c} g^{b d} \nabla_d \omega \tag{2}
$$
Recall that $\omega = \log \Omega$ and introduce the notation:
$$
\Upsilon_a := \nabla_a \log \Omega \tag{3}
$$
It is also customary in Riemannian geometry to lower and raise the indices using the Riemannian metric $g_{a b}$ and its inverse $g^{a b}$ without mentioning this explicitly, so, for instance, $\Upsilon^a = g^{a b} \Upsilon_b$.
Keeping in mind that $X^a$ is arbitrary, we can factor it away from equation $(2)$, and renaming along the way $Y$ to $X$, we obtain:
$$
\widehat{\nabla}_a X^b = \nabla_a X^b + \Upsilon_a X^b + \delta_a{}^b \Upsilon_c X^c - X_a \Upsilon^b \tag{4}
$$
Here $\delta^a{}_b$ is the identity operator, also known as the Kronecker symbol.
It is convenient to introduce the following difference tensor
$$
S_a{}^b{}_c := \Upsilon_a \delta^b{}_c + \delta_a{}^b \Upsilon_c - \Upsilon^b g_{a c} \tag{5}
$$
Using $(5)$ we can rewrite $(4)$ simply as
$$
\widehat{\nabla}_a X^b = \nabla_a X^b  + S_a{}^b{}_c X^c \tag{6}
$$
For any tensor $A^{a b}$, we can calculate the conformal rescaling rule of its covariant derivative using the Leibniz rule in disguise:
$$
\widehat{\nabla}_a A^{b c} = \nabla_a A^{b c} + S_a{}^b{}_d A^{d c} + S_a{}^c{}_d A^{b d} \tag{7}
$$
The divergence $\nabla_a A^{a b}$ is then rescaled according to the identity
$$
\widehat{\nabla}_a A^{a b} = \nabla_a A^{a b} + S_a{}^a{}_c A^{c b} + S_a{}^b{}_c A^{a c} \tag{8}
$$
We can simplify the terms of the above equation, using $(5)$, as follows:
$$
S_a{}^a{}_c A^{c b} = (\Upsilon_c + n \Upsilon_c - \Upsilon_c) A^{c b} = n \Upsilon_c A^{c b} \tag{9}
$$
and
$$
S_a{}^b{}_c A^{a c} = (\Upsilon_a \delta^b{}_c + \delta_a{}^b \Upsilon_c - \Upsilon^b g_{a c}) A^{a c} = \Upsilon_a A^{a b} + \Upsilon_c A^{b c} - \Upsilon^b \mathrm{tr} A \tag{10}
$$
And now, if $A^{a b}$ is symmetric ($A^{a b} = A^{b a}$) and trace-free ($\mathrm{tr}A = 0$), the last equation simplifies even further
$$
S_a{}^b{}_c A^{a c} = 2 \Upsilon_c A^{c b} \tag{11}
$$
Thus, for a symmetric trace-free $A^{a b}$, we have the following formula for the conformal rescaling of its divergence:
$$
\widehat{\nabla}_a A^{a b} = \nabla_a A^{a b} + (n + 2) \Upsilon_a A^{a b}  \tag{12} 
$$
This means that for any (trace-free symmetric) tensor $A^{a b}$ its divergence is not conformally invariant, unless $A^{a b}$ is identically zero (everywhere vanishing).
We can improve our situation by considering weighted versions of $A^{a b}$. Namely, let $w \in \mathbb{R}$ be any real number, so that we can calculate
$$
\widehat{\nabla}_a (\Omega^w A^{a b})
= (\widehat{\nabla}_a \Omega^w) A^{a b}  + \Omega^w \widehat{\nabla}_a A^{a b} \\
= (w \Omega^{w - 1} \nabla_a \Omega) A^{a b} + \Omega^w (\nabla_a A^{a b} + (n + 2) \Upsilon_a A^{a b}) \\
= (w \Omega^{w} \nabla_a \log \Omega) A^{a b} + \Omega^w (\nabla_a A^{a b} + (n + 2) \Upsilon_a A^{a b}) \\
= \Omega^w \big( \nabla_a A^{a b} + (n + w + 2) \Upsilon_a A^{a b} \big) \tag{13}
$$
The weight is something which we can choose, so that if $\dim M = n$, choosing $w = - n - 2$, we can achieve the following conformal covariance property:
$$
\widehat{\nabla}_a (\Omega^{ - n - 2} A^{a b}) = \Omega^{ - n - 2} \nabla_a A^{a b} \tag{14}
$$
In particular, if $n = 3$, we can use the weight $-5$, and $\psi^{- 2} = \Omega$, to recover the result in OP.
