Consider the an infinitesimal coordinate transformation $x^{\mu}\mapsto x'^{\mu}:=x^{\mu}+\epsilon^{\mu}(x)$. The appropriate changes in the metric manifest as: $$g'_{\mu\nu}=\frac{\partial x^{\alpha}}{\partial x'^{\mu}}\frac{\partial x^{\beta}}{\partial x'^{\nu}}g_{\alpha\beta}= g_{\mu\nu}-(\partial_{\mu}\epsilon_{\nu}+\partial_{\nu}\epsilon_{\mu})+\mathcal{O}(\epsilon^2).$$ Now, for the transformation above to qualify as a conformal transformation, the $\mathcal{O}(\epsilon^1)$-term must be proportional to the metric $g_{\mu\nu}$ and the proportionality factor can be determined by taking a trace of the resultant expression. All in all, $$\partial_{\mu}\epsilon_{\nu}+\partial_{\nu}\epsilon_{\mu}\overset{!}{=}\frac{2}{D}\partial\cdot\epsilon g_{\mu\nu},\quad(*)$$ where $D:=n+m$ is the dimensionality of the semi-Riemannian manifold $\mathbb{R}^{n,m}$.
The main problem I am facing here is using (*) to get to, $$(g_{\mu\nu}\square +(D-2)\partial_{\mu}\partial_{\nu})(\partial\cdot\epsilon)=0,\quad(**)$$ where $\square$ is the d'Alembertian [1]. According to a lecture by Prof. Tobias Osborne, "it is due to the equality of mixed partial derivatives [...] up to third order ..." [2].
The method I have adopted is to evaluate $\partial_{\alpha}\partial_{\beta}\times(*)$ and note that the RHS(*) is $\propto g_{\mu\nu}$ and therefore it is sufficient to consider $\partial_{\mu}\epsilon_{\nu}=\partial\cdot\epsilon g_{\mu\nu}$ to get $$(2\partial_{\alpha}\partial_{\beta}\partial\cdot\epsilon)g_{\mu\nu}=\frac{2}{D}\partial_{\alpha}\partial_{\beta}\partial\cdot\epsilon g_{\mu\nu},$$ which appears a fairly bogus result. How do I get to (**)?
[1] https://arxiv.org/abs/hep-th/9108028
[2] https://youtu.be/NGYX6gtObec?list=PLDfPUNusx1Ep5g1jIKXqpNX_t_Zz-kAlQ&t=2685
