Let $(M,g)$ be a semi-Riemannian manifold. A Conformal Killing Vector (CKV) $Y$ is any vector field obeying the Conformal Killing Equation: $$\nabla_\mu Y_\nu+\nabla_\nu Y_\mu=g_{\mu\nu}\dfrac{2}{d} \nabla\cdot Y,\quad \nabla\cdot Y\equiv \nabla_\mu Y^\mu\tag{1}.$$
Contracting the equation with $\nabla^\mu\nabla^\nu$ and using the Riemann tensor I found this equation can be written as $$\left(2-\frac{2}{d}\right)\Box(\nabla\cdot Y)+\nabla^\mu (R^{\nu}_{\phantom{\nu}\rho\nu\mu}Y^\rho)+\nabla^\nu (R^\mu_{\phantom{\mu}\rho\nu\mu}Y^\rho)+R^{\nu}_{\phantom \nu \rho \mu\nu}\nabla^\rho Y^\mu+R^{\mu}_{\phantom \mu \rho \mu\nu}\nabla^\nu Y^\rho=0,\quad \Box \equiv \nabla^\mu \nabla_\mu\tag{2}$$
In maximally symmetric spaces, using the form of $R^\sigma_{\phantom{\sigma}\rho\mu\nu}$, I managed to show that $$\Box(\nabla\cdot Y)=0.\tag{3}$$
This seems to be a very general property of CKVs.
Now let us consider the sphere $(S^2,\gamma)$ with round metric $\gamma$. We shall denote the indices with uppercase latin letters and let $D$ be the Levi-Civita connection. Also I write $D^2$ for $\Box$. In the analysis of the $\mathfrak{bms}_4$ algebra there appears one such CKV $Y$ on $S^2$. In trying to reproduce the results from the literature I came to a point in which to obtain the known results it would be necessary to employ $$D^2(D\cdot Y)=-2D\cdot Y\tag{4}.$$
Now it would seem that this is impossible in view of (3). If $Y$ is a CKV the LHS of (4) must be zero since $(S^2,\gamma)$ is maximally symmetric!
It turns out, however, that there is mention to (4) in the literature. For example in the solutions section of "Lectures on the IR Structure of Gravity and Gauge Theory", in page (132) it is said right bellow Eq. (9.0.180)
Note that $$\gamma^{AB}\delta C_{AB} = \frac{u}{2}(-4D\cdot Y-2D^2D\cdot Y)+C^{AB}(D_AY_B+D_BY_A)=0,\tag{9.0.181}$$ where $\pmb{D^2D\cdot Y=-2D\cdot Y}$, and $C^{AB}(D_AY_B+D_BY_A)=C^{AB}\gamma_{AB}D\cdot Y$. Hence, $C_{AB}$ remains traceless under the diffeomorphism generated by $\zeta$.
So now I'm quite confused because this is saying that (4) is true. But since $Y$ is a CKV (3) should be true and both are incompatible unless $D\cdot Y=0$ which is certainly not the case in general.
What am I missing here? I am quite sure I'm missing something very basic here. Is this some issue with notation, so that the LHS of (4) is not the same as the LHS of (3)? What is the actual relation between (3) and (4)?
Your Equation (2), and hence Equation (3), is incorrect. The correct version of Equation (2) is $$ \left( 2 - \frac{2}{d}\right)\Box(\nabla\cdot Y)+ 2\nabla^\mu\left(R^\nu{}_{\rho\nu\mu}Y^\rho\right) + R^\nu{}_{\rho\mu\nu}\nabla^\rho Y^\mu + R^\mu{}_{\rho\mu\nu}\nabla^\nu Y^\rho = 0 . $$ For what it's worth, this immediately simplifies to $$ \left(2-\frac{2}{d}\right)\Box(\nabla\cdot Y) + 2\nabla^\mu(R_{\rho\mu}Y^\rho) $$ on all pseudo-Riemannian manifolds, which is a much nicer formula to look at. Evaluating on the two-sphere is now consistent with your Equation (4).