Introduction
A conformal Killing vector field $\xi^\mu$ is one that satisfies $$\mathcal{L}_\xi g_{\mu\nu}\equiv\nabla_\mu\xi_\nu+\nabla_\nu\xi_\mu=\kappa g_{\mu\nu}$$ for some scalar field $\kappa$. Clearly we have $\kappa=(2/n)\nabla_\mu\xi^\mu$ where $n$ is the dimensionality of the manifold.
I am starting to read Schottenloher's text book on Conformal Field Theory and he mainly focuses on flat manifolds. In particular a very useful theorem that he proves for flat manifolds is that the conformal $factor$, $\kappa$, satisfies $$(n-2)\partial_\alpha\partial_\beta \kappa + g_{\alpha\beta}\partial^\mu\partial_\mu\kappa=0$$
He also says that this generalises to curved manifolds as $$(n-2)\nabla_\alpha\nabla_\beta \kappa + g_{\alpha\beta}\nabla^\mu\nabla_\mu\kappa=0$$
THE QUESTION
How do you prove the generalised condition for curved space times? $$(n-2)\nabla_\alpha\nabla_\beta \kappa + g_{\alpha\beta}\nabla^\mu\nabla_\mu\kappa=0$$
MY PROGRESS
Inspired by Schottenloher's proof for flat spacetimes I write $$(n-2)\nabla_\mu\nabla_\nu \kappa + g_{\mu\nu}\nabla^\alpha\nabla_\alpha\kappa=g^{\alpha\beta}\Big(g_{\mu\nu}\nabla_\alpha\nabla_\beta+g_{\alpha\beta}\nabla_\mu\nabla_\nu-g_{\alpha\mu}\nabla_\beta\nabla_\nu-g_{\beta\nu}\nabla_\mu\nabla_\alpha\Big)\kappa$$ $$=g^{\alpha\beta}\Big[\big(\nabla_\alpha \nabla_\beta \nabla_\nu - \nabla_\beta \nabla_\nu \nabla_\alpha\big)\xi_\mu + \big(\nabla_\alpha \nabla_\beta \nabla_\mu - \nabla_\mu \nabla_\alpha \nabla_\alpha\big)\xi_\nu$$ $$\hspace{12mm} + \big(\nabla_\mu \nabla_\nu \nabla_\beta - \nabla_\beta \nabla_\nu \nabla_\mu\big)\xi_\alpha + \big(\nabla_\mu \nabla_\nu \nabla_\alpha - \nabla_\mu \nabla_\alpha \nabla_\nu\big)\xi_\beta\Big]$$ $$=g^{\alpha\beta}\Big[\nabla_\beta\big(R_{\alpha\nu\mu\sigma}\xi^\sigma\big) + \nabla_\alpha\big(R_{\beta\mu\nu\sigma}\xi^\sigma\big)$$ $$+\nabla_\beta\big(R_{\mu\nu\alpha\sigma}\xi^\sigma\big)+\nabla_\mu\big(R_{\nu\beta\alpha\sigma}\xi^\sigma\big)+\nabla_\mu\big(R_{\nu\alpha\beta\sigma}\xi^\sigma\big)$$ $$+R_{\alpha\beta\nu\sigma}\nabla^\sigma\xi_\mu+R_{\alpha\beta\mu\sigma}\nabla_\nu\xi^\sigma+R_{\alpha\mu\beta\sigma}\nabla^\sigma\xi_\nu$$ $$+R_{\alpha\mu\nu\sigma}\nabla_\beta\xi^\sigma+R_{\mu\beta\nu\sigma}\nabla^\sigma\xi_\alpha+R_{\mu\beta\alpha\sigma}\nabla_\nu\xi^\sigma\Big]$$ $$=2\nabla_\alpha\big(R^\alpha_{\hspace{2mm}\nu\mu\sigma}\xi^\sigma\big)-2\nabla_\mu\big(R_{\nu\alpha}\xi^\alpha\big)+R_{\mu\alpha}\big(\nabla^\alpha\xi_\nu-\nabla_\nu\xi^\alpha\big)+R^\alpha_{\hspace{2mm}\mu\nu\sigma}\big(\nabla_\alpha\xi^\sigma-\nabla^\sigma\xi_\alpha\big)$$
I know I should use the Bianchi identity and a whole lot of Riemann's symmetry properties, but how?
Also, I didn't find any source that has this result in it, if you introduce some, I would also consider it as an answer; Thanks!