I've always been very uncomfortable with the following two conformal Ward identities [c.f. Di Francesco, pg 107, Eqs. (4.66) and (4.67)]: $$ \langle (T^{\rho \nu} - T^{\nu \rho})X \rangle = -i \sum_i \delta(\boldsymbol{x} - \boldsymbol{x}_i) S^{\nu \rho}_i \langle X \rangle, \quad \langle T^{\mu}_{\mu} X \rangle = -\sum_i \delta(\boldsymbol{x} - \boldsymbol{x}_i) \Delta_i \langle X \rangle $$ where $X = \phi_1(\boldsymbol{x}_1) \ldots \phi_n(\boldsymbol{x}_n)$ is a product of local fields with spin generators $S^{\nu \rho}_i$ and scaling dimensions $\Delta_i$. In particular, I've always been confused how the identities $T^{\mu \nu} = T^{\nu \mu}$ and $T^{\mu}_{\mu} = 0$ can hold as operator equations, and yet be violated at insertions of other fields. (In the more familiar Ward identities for a conserved current, $\partial_{\mu} \langle j^{\mu} X \rangle = \ldots$, I am much more comfortable: $\partial_0$ should be interpreted as acting on the time-ordering symbol, so $\partial_{\mu} j^{\mu} = 0$ holding as an operator equation isn't inconsistent with $\partial_{\mu} \langle j^{\mu} X \rangle$ being nonvanishing at field insertions.)
To try and alleviate my confusion, I want to verify these equations in the simple case of a 2d free scalar field $\phi$. Here is my attempt. First of all, in 2d complex coordinates, the above two equations can be rewritten as [c.f. Di Francesco, pg 120, Eqs. (5.37) and (5.38)] $$ 2\langle T_{\bar{z} z} X \rangle = -\sum_i \delta(\boldsymbol{x} - \boldsymbol{x}_i) h_i \langle X \rangle, \quad 2\langle T_{z\bar{z}} X \rangle = -\sum_i \delta(\boldsymbol{x} - \boldsymbol{x}_i) \bar{h}_i \langle X \rangle $$ where $h_i = \frac{1}{2}(\Delta_i + s_i)$ and $\bar{h}_i = \frac{1}{2}(\Delta_i - s_i)$ are the conformal dimensions of $\phi_i$. In the case of a scalar field, the stress tensor is $T_{\mu \nu} = \partial_{\mu} \phi \partial_{\nu} \phi - \frac{1}{2} \eta_{\mu \nu} (\partial \phi)^2$, so the off-diagonal components of $T_{\mu \nu}$ in complex coordinates are $$ T_{z \bar{z}} = \partial \phi \bar{\partial} \phi - \frac{1}{2} (\partial \phi \bar{\partial}\phi + \bar{\partial}\phi \partial \phi) = \frac{1}{2} (\partial \phi \bar{\partial}\phi - \bar{\partial} \phi \partial \phi) $$ and an identical calculation seems to give $T_{\bar{z} z} = -T_{z \bar{z}}$. It's very tempting to declare that $T_{z \bar{z}} = T_{\bar{z} z} = 0$, but this is clearly incorrect; I believe the issue has to do with the necessity of regularizing $T_{\mu \nu}$, so let me assume that the above expression is appropriately defined by point-splitting.
Now, to verify a particular case of the above Ward identities, let's take $X = \partial \phi(w_1) \partial \phi(w_2)$, so that $h_1 = h_2 = 1$ and $\bar{h}_1 = \bar{h}_2 = 0$. The result we want to get is $$ \langle 2 T_{\bar{z} z}(\boldsymbol{y}) X \rangle = \left[ \delta(\boldsymbol{y} - \boldsymbol{x}_1) + \delta(\boldsymbol{y} - \boldsymbol{x}_2) \right] \frac{1}{4\pi (w_1 - w_2)^2}, \quad \langle 2 T_{z \bar{z}}(\boldsymbol{x}) X \rangle = 0 \tag{*} $$ where in the above, $\boldsymbol{x}_1 = (x^0_1, x^1_1)$ and $w_1 = x^0_1 + i x^1_1$ (and similarly for $\boldsymbol{x}_2$ and $w_2$). To get the same result from the lefthand side, we need to evaluate correlation functions of the form (ignoring constant prefactors for the moment) $$ \langle T_{\bar{z} z}(\boldsymbol{y}) \partial \phi(w_1) \partial \phi(w_2) \rangle \sim \lim_{z_1, z_2 \to z} \left[ \frac{\delta(\boldsymbol{y}_2 - \boldsymbol{x}_2)}{(z_1 - w_1)^2} + \frac{\delta(\boldsymbol{y}_2 - \boldsymbol{x}_1)}{(z_1 - w_2)^2} - \frac{\delta(\boldsymbol{y}_1 - \boldsymbol{x}_2)}{(z_2 - w_1)^2} - \frac{\delta(\boldsymbol{y}_1 - \boldsymbol{x}_1)}{(z_2 - w_2)^2} \right] $$ where I've used $\langle \partial \phi(z) \bar{\partial} \phi(\bar{w}) \rangle = \frac{1}{4} \delta(\boldsymbol{y} - \boldsymbol{x})$, and $z$ is the complex coordinate for $\boldsymbol{y}$. At this stage, I already feel as though something has gone horribly wrong: if I performed an analogous computation with $T_{z \bar{z}}$ instead of $T_{\bar{z} z}$, I would get exactly minus this result. So, whether this expression is zero or nonzero upon taking the limit $z_1, z_2 \to z$, I get the incorrect result for one of the two Ward identities!
So, here is the question: what is the correct way to evaluate the lefthand side of the two equations $(*)$, in order to get the desired result from the Ward identities?