For local observables $\{O_i(x_i)\}^n_{i = 1}$, one defines the Ward identity as
$$\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O(x_i)\rangle = \sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O(x_n)\rangle\tag{1}$$
Under the local variation of the field $$\Psi(x)\mapsto \Psi(x)+\epsilon\delta\Psi(x)\tag{2}$$
Whenever $\delta_{\epsilon}\mathcal{D}[\Psi(x)] = 0$, $\delta_{\epsilon} Z[\Psi(x)] = 0$, one obtains the following variation
$$\begin{align}\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle &= \delta_{\epsilon}\biggr(\frac{\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}}{Z[\Psi(x)]}\biggr)\\&= \frac{\delta_{\epsilon}\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)Z[\Psi(x)] – \delta_{\epsilon} Z[\Psi(x)]\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)}{Z^2[\Psi(x)]}\\&= \frac{1}{Z[\Psi(x)]}\delta_{\epsilon}\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)-\underbrace{\frac{1}{Z[\Psi(x)]}\delta_{\epsilon} Z[\Psi(x)]\langle \prod^n_{i = 1}O_i(x_i)\rangle}_{ = 0} \\ &= \sum^n_{i = 1}\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\langle\delta_{\epsilon} S[\Psi(x)]\prod^n_{i = 1}O_i(x_i)\rangle\\&= \sum^n_{i = 1}\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\int_{M} d^n x\epsilon(x)\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle\end{align}\tag{3}$$
Where we assumed that $$Z[\Psi(x)] = \displaystyle\int D[\Psi(x)]e^{\frac{i}{\hbar}S[\Psi(x)]}\tag{4}$$
Further recalling that $$\int_{V}O_i(x)\delta(x-x_i)d^nx = \begin{cases}O_i(x_i), x_i\in V\\ 0, x_i\notin V\end{cases}\tag{5}$$
We observe that under arbitrary spacetime-dependent parameter $\epsilon : M\rightarrow \mathbb{R}$,
$$\begin{align}\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle &= \int_{M}d^nx\sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\int_{M} d^nx\epsilon(x) \partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle\end{align}\tag{6}$$
From which we conclude that whenever $\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle = 0$,
$$-i\hbar\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O(x_i)\rangle = \sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O(x_n)\rangle\tag{7}$$
However, I am not sure how the Dirac delta function $\delta(x-x_i)$ arises in the Ward identity in $(1)$. Can you provide clarification?
Best Answer
It might be completely wrong, if it sounds fishy, please tell me.
You have to keep in mind the difference between $\delta\langle O\rangle$ and $\langle \delta O\rangle$. The former is 0 by definition of the invariance of the system, while the latter is not (somehow it accounts for the variation of the observable that compensate the variation of the probability density so that the total variation is 0). This implies that your equation (2) (which corresponds to the former case) is equal to 0. Moreover, the variation here: $$O_i(x_i)\mapsto O_i(x_i)+\epsilon\delta O_i(x_\color{red}i\color{black})$$ is with respect to the space time variable you're considering (and not an arbitrary $x$). These two things lead us to: $$\sum_{i = 1}\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O_n(x_n)\rangle=-i\int_{M} d^n x\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle$$
This gives us, using the delta identity: $$\int_M d^nx\sum_{i = 1}\delta(x_i - x)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O_n(x_n)\rangle=-i\int_{M} d^n x\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle$$.
Since this integral holds for every variation, the equality holds for the objects inside the integrals:
$$\sum_{i = 1}\delta(x_i - x)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O_n(x_n)\rangle=-i\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle$$
You can make the last argument clearer by clearly using a $x$ dependent infinitesimal: $\epsilon(x)$ (which we omitted here). The argument is the same for the derivation of the Euler-Lagrange equations.
Edit: Concerning your question in the comments. I don't know a damn thing about QFT :) ! So my answer would be to ask directly with another question (I'm also interested by this). My (classical mechanics) feeling would be that, your infinitesimal spacetime variations leaving unchanged your average values belong to a Lie algebra and thus, through the exponential map can form the full Lie group associated to this algebra. This means that, even if you only consider infinitesimal variation, you can reconstruct, FROM THESE INFINITESIMAL variations, any finite element of the Lie group (through exponentiation). What this means, to me, is that there is as much information in the $\epsilon$ Ward identity that there would be in an $\epsilon^2$ (or even finite variation) Ward identity.
Here are a PSE link about these ideas of Lie algebras. Here is also a PSE post, with an answer by ACuriousMind who says at the end:
Here is an Article saying that my argument does not always hold (albeit, in a rather specific context):