Conformal Ward Identity

I'm trying to prove: (Exercise from "TASI Lectures on the Conformal Bootstrap" by David Simmons-Duffin)

$$\partial^{\mu} \langle T^{\mu \nu} O_{1}(x_{1}) \dots O(x_{n}) \rangle = \sum_{i = 1}^{N}\delta(x - x_{i}) \partial^{\nu} \langle O_{1}(x_{1}) \dots O(x_{n}) \rangle$$

The theory has to be invariant under a diffeomorphism of $$S[g, \phi]$$, thus we have the following equalities:

1. $$x^{\prime\mu} = x^{\mu} + \epsilon^{\mu}(x)$$
2. $$\dfrac{\partial x^{\prime\mu}}{\partial x^{\nu}} = \delta^{\mu}_{\; \nu} + \partial_{\nu}\epsilon^{\mu}$$
3. $$g^{\prime\mu \nu} = \dfrac{\partial x^{\prime\mu}}{\partial x^{\sigma}}\dfrac{\partial x^{\prime\nu}}{\partial x^{\sigma}}g^{\sigma \sigma} = g^{\mu \nu} + \left( \partial^{\mu}\epsilon^{\nu} + \partial^{\nu}\epsilon^{\mu}\right) + O \left(\epsilon^2 \right)$$

As this is a symmetry, I wrote:

$$S[\phi^{\prime}] = S[\phi], \quad D\phi^{\prime} = D\phi, \quad \langle O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n}) \rangle_{g} = \langle O_{1}(x_{1}) \dots O(x_{n}) \rangle_{g}$$

Doing the transformation inside the path integral we modify $$O(x_{1}^{\prime}) = O(x_{1}^{\prime} + \epsilon) = O(x_{1}) +\epsilon^{\nu}\partial_{\mu}O(x_{1})$$

So:

$$\langle O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n}) \rangle_{g} = \int D\phi^{\prime} O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n}) e^{-S[g^{\prime}, \phi^{\prime}]} \\ =\int D\phi O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n}) e^{-S[g, \phi] + \delta g^{\mu \nu} \frac{\delta S[g, \phi] }{\delta g^{\mu \nu}}}$$ Expanding the exponential

$$=\int D\phi O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n})\left[ 1 - \delta g^{\mu \nu} \frac{\delta S[g, \phi] }{\delta g^{\mu \nu}} \right] e^{-S[g, \phi]} \\ = \int D\phi \left[O(x_{1}) +\epsilon^{\nu}\partial_{\mu}O(x_{1}) \right] \dots \left[O(x_{n}) +\epsilon^{\nu}\partial_{\mu}O(x_{n}) \right]\left[ 1 - \delta g^{\mu \nu} \frac{\delta S[g, \phi] }{\delta g^{\mu \nu}} \right] e^{-S[g, \phi]}$$ To mantain $$\langle O^{\prime}_{1}(x_{1}) \dots O^{\prime}(x_{n}) \rangle_{g} = \langle O_{1}(x_{1}) \dots O(x_{n}) \rangle_{g}$$ and neglecting higher orders, we get the following condition:

$$\langle \frac{\delta S[g, \phi] }{\delta g^{\mu \nu}} O_{1}(x_{1}) \dots O(x_{n}) \rangle = - \sum_{i = 1}^{N} \epsilon^{\nu} \partial_{\mu} \langle O_{1}(x_{1}) \dots O(x_{n}) \rangle$$

That is not the conformal Ward identity. Am I missing something?

Thank you very much for your help!

In my opinion, it is more instructive to derive Ward-Takahashi identitities for the general case, $$\partial_{\mu} \left\langle J_a^{\mu}(x) \Phi(x_1) \ldots \Phi(x_n)\right\rangle = -\sum_i \delta(x - x_i) \left\langle \Phi(x_1) \ldots X_a(\Phi)(x_i) \ldots \Phi(x_n)\right\rangle,$$ with $$\Phi(x) \to \Phi(x) + \epsilon_a(x) X_a(\Phi)(x)$$, and then apply it to specific symmetries, but okay.
First of all, why would $$O \to O + \epsilon^{\nu} \partial_{\mu} O$$? Let's forget for a second that that's essentially a Taylor expansion. The fact that on the right-hand side you have a spacetime scalar $$+$$ a spacetime tensor should have alarmed you that something isn't right. Of course, you should have written $$\epsilon_{\mu}\partial^{\mu} O$$ instead. Furthermore, notice that in each term the argument will be different. So the right-hand side in your last equation shall read $$-\sum_i \epsilon_{\mu}(x_i) \frac{\partial}{\partial x_i^{\mu}} \left\langle O(x_1) \ldots O(x_n) \right\rangle.$$ On the left-hand side, you forgot to multiply by $$\delta g^{\mu\nu}$$ (as well as the minus sign?). But let's actually not bother with $$\delta g$$ at all and simply write the change of action as $$\delta S$$. Putting everything together we then obtain $$\left\langle \delta S\, O(x_1) \ldots O(x_n) \right\rangle = \sum_i \epsilon_{\mu}(x_i) \frac{\partial}{\partial x_i^{\mu}} \left\langle O(x_1) \ldots O(x_n) \right\rangle.$$ Now, you probably remember that for symmetry transformations the change of action takes the form $$\delta S = -\int_x \epsilon_k \partial_{\mu} J_k^{\mu} = \int J_k^{\mu} \partial_{\mu} \epsilon_k$$, with $$J_k^{\mu}$$ being conserved currents. For translation symmetry, one has a current for each of the translations $$\epsilon_{\nu}$$, with $$\nu = 0, 1, 2, 3$$, and the change in action reads $$\delta S = \int_x T^{\mu\nu}(x) \partial_{\mu}\epsilon_{\nu}(x)$$ (derive this result if you haven't done this before!). Thus, $$\int_x \partial_{\mu}\epsilon_{\nu}(x) \left\langle T^{\mu\nu}(x) O(x_1) \ldots O(x_n) \right\rangle = \int_x \sum_i \delta(x - x_i) \epsilon_{\nu}(x) \frac{\partial}{\partial x^{\nu}} \left\langle O(x_1) \ldots O(x_n) \right\rangle,$$ where we have made some obvious algebraic manipulations on the right-hand side. Finally, upon integrating the left-hand side by parts, $$-\int_x \epsilon_{\nu}(x) \partial_{\mu} \left\langle T^{\mu\nu}(x) O(x_1) \ldots O(x_n) \right\rangle = \int_x \epsilon_{\nu}(x) \sum_i \delta(x - x_i) \frac{\partial}{\partial x^{\nu}} \left\langle O(x_1) \ldots O(x_n) \right\rangle,$$ and noting that the above identity should hold for any infinitesimal $$\epsilon_{\mu}$$ we conclude that we may write the following local relation: $$\partial_{\mu} \left\langle T^{\mu\nu}(x) O(x_1) \ldots O(x_n) \right\rangle = -\sum_i \delta(x - x_i) \frac{\partial}{\partial x^{\nu}} \left\langle O(x_1) \ldots O(x_n) \right\rangle,$$ which is exactly the desired Ward identity corresponding to translation symmetry (modulo sign errors and such).