Let $(M, g)$ be a compact Riemannian manifold of dimension $n \geq 4$.

In Besse's book, the Bach tensor is defined as the gradient of the functional

$$ SW(g) = \int_M |W(g)|_g^2 d\mu^g $$

meaning that, for any symmetric 2-tensor field,

$$ SW(g + t h) = SW(g) + t \int_M \langle B, h\rangle_g d\mu^g + O(t^2) $$

As a consequence, this symmetric 2-tensor is

- Divergence-free as one has $SW(\phi_t^* g) = SW(g)$ for any 1-parameter family of diffeomorphisms $\phi_t$. So choosing for $\phi_t$ the flow of some vector field $X$, we get $\phi_t^* g = g + t L_X g + O(t^2)$. As $$ SW(\phi_t^* g) = SW(g) + t \int_M \langle B, L_X g\rangle_g d\mu^g + O(t^2), $$ we conclude that $$ \int_M \langle \mathrm{div} B, X\rangle_g d\mu^g = \int_M \langle B, L_X g\rangle_g d\mu^g = 0 $$ for any vector field $X$.
- Trace-free in dimension 4.This comes from the fact that $SW$ is conformally invariant in dimension $4$ so one can repeat the previous argument with the 1-parameter family of metric $g_t = e^{\lambda u} g$ (for an arbitrary function $u$).

Now the formula given for the Bach tensor (see e.g. the Wikipedia page of the Bach tensor, but this is very similar to the formula in Besse's book) is

$$ B_{ab} = P^{cd} W_{acbd} + \nabla^c (\nabla_c P_{ab} - \nabla_a P_{cb}) $$ (here $P$ is the Schouten tensor)

This formula for $B$ gives a trace-free tensor in any dimension. This follows immediately as the Weyl tensor is trace-free and as the second term is the divergence of the Cotton tensor, so, up to some coefficient, this is the double divergence of the Weyl tensor. This is independent of the dimension.

The divergence, however, can be computed explicitely as $$ \nabla^a B_{ab} = P^{cd} (\nabla_c P_{db} - \nabla_b P_{cd}) $$ which is a priori non-zero.

I do not understand where my mistake is. Any help would be appreciated.