Suppose we have $\theta$-field in QCD (in a special case of constant $\theta$ it reduces to ordinary $\theta$ parameter): $$ \tag 1 Z_{\theta} = \int D[\psi_{QCD}]e^{iS}, $$ where $$ \tag 2 S = \int d^4x \left( \bar{q}\gamma_{\mu}D^{\mu}q - \bar{q}^{i}M_{ij}q^{j} - \frac{1}{4g_{s}^2}GG + \theta G\tilde{G} + \text{gauge}+\text{ghost} \right), $$ and $\psi_{QCD}$ are all fields which take part in QCD interactions.
I need to calculate $\theta$ parameter vacuum energy, more precisely, the quadratic term. In principle, from $(1), (2)$ it can be evaluated by computing quantity $$ \tag 3 \frac{1}{2V_{4}}\int d^{4}xd^{4}y \left(\frac{\delta^{2}Z}{\delta \theta (x)\delta \theta (y)}\right)_{\theta = 0} \equiv \frac{1}{2V_{4}}\int d^4x d^4y \langle 0| G\tilde{G}(x)G\tilde{G}(y)|0\rangle_{\theta = 0} \equiv \frac{1}{2}\kappa (0), $$ where $\kappa (0)$ is QCDs topological susceptibility. From $(3)$ it seems that coordinate dependence of $\theta$ doesn't play role; in particular, energy of constant $\theta$ term at zero momentum will be the same as of variable one.
But I've made following trick for evaluating $(3)$: I've made local chiral transformation $q \to e^{iC\gamma_{5}\theta (x)}q$, where $C$ is constant, and have get from $(1), (2)$ generating functional in a form $$ Z_{\theta}[...] = \int D[\psi_{QCD},\theta]e^{iS{'}}, $$ where $$ \tag 5 S{'} = \int d^4x\left( \bar{q}\gamma_{\mu}D^{\mu}q - \bar{q}^{i}M_{ij}e^{2iC\gamma_{5}\theta(x)}q^{j} - \frac{1}{4g_{s}^2}GG + Ci\partial_{\mu}\theta J^{\mu}_{5}+ \text{gauge}+\text{ghost} \right) $$ In the case of constant $\theta$ the fourth term is absent. Now I can substitute $(5)$ into $(3)$. After tedious calculations I will get that in compare to the case of constant $\theta$ I will also get extra term $$ E_{2} = \frac{1}{2V_{4}}\times 8C^2\int d^4x d^4y\delta (x - y) \langle 0|\bar{q}^{i}M_{ij}q^{j} |0\rangle , $$ which comes from the initial correlator $$ \frac{1}{2V_{4}}4C^2 \int d^4y d^4x \langle 0|T\left( \partial_{\mu}J^{\mu}_{5}(x)\bar{q}^{i}\gamma_{5}M_{ij}q^{j}(y)\right)| 0\rangle $$ and almost coincides with corresponding term $E_{1}$ for constant $\theta$: $$ E = E_{1} + E_{2} $$ So or I've made the mistake in the statement that the vacuum energy of $\theta$ field at zero momentum doesn't depend on its coordinate dependence, or I've made incorrect calculations, or my result is true.
Can anyone help me with this problem?
An edit
It seems that calculations given above contain serious error. It is closely related to this question, Time-ordered product vs path integral, because in the calculations above I first used path integral approach, while later I used canonical operator formalism. The problem is that $(3)$ is true in path integral approach but isn't true in canonical operator formalism, while the next expressions are true in canonical operator formalism, but aren't true in path integral approach. I need to use path integral formalism if I want to compute explicit form of $(3)$; in this approach I may simply throw out derivatives, and contact terms don't arise; otherwise for working in operator formalism I need to derive the correct form of $(3)$ which includes contact (so-called Schwinger) terms, and they will coincide boundary terms like $(2)$.
