Why is there no gauge-invariant local operator in GR? I have a hard time understanding why the bulk locality is a question. I know some operator which depends on a particular coordinate $x$, $O(x)$, and its correlation function like $ \langle O(x)O(y) \rangle$ for example can’t be diffeomorphism invariant as the operator itself becomes different in some other coordinate patch, while nonlocal operators (i.e. operators integrated over whole spacetime region) is gauge invariant.
Then, why don’t we think an operator $O(p)$ where $p \in M$ on a spacetime manifold $M$? At the manifold level, we are yet to introduce a coordinate system and thus this kind of operator should be naturally gauge invariant.
 A: For gauge theories, observables (schematically) are functions of the form
\begin{equation}
\mbox{observable}\colon
\frac{\mbox{phase space}}{\mbox{gauge symmetries}}\to\mathbb{R},
\end{equation}
which maps distinct physical states ["gauge orbits"] to "some number"
(see arXiv:hep-th/0203027 for a review of the geometry of gauge systems, or the first 5 chapters of Henneaux and Teitelboim's Quantization of Gauge Systems). This condition is satisfied if the observable commutes with the Poisson bracket for each generator of the gauge symmetry, i.e., every first-class constraint. General relativity has the diffeomorphism group for its gauge symmetry, which causes complications: the algebra of constraints is infinite-dimensional.
Torre (arXiv:gr-qc/9306030) provides a "no-go" theorem denying the existence of  any nontrivial local and diffeomorphism-invariant observables. This is on "closed" universes, i.e., the spatial hypersurfaces decomposition of the spacetime into $\mathbb{R}\times\Sigma$ has the space-like hypersurface $\Sigma$ be compact without boundary.
We can try to construct a class of observables by considering integrals over spacetime of scalar functions, but this family of observables do not posses a local interpretation. The scalar must commute with the generators of the diffeomorphism constraints. For a given scalar field $\phi(x)$, the diffeomorphism group acts as $\partial_{\mu}\phi(x)$ which vanishes if and only if $\phi(x)=\phi_{0}$ is a constant. Torre has shown general relativistic observables must include an infinite number of derivatives and hence are very nonlocal.
The proposed routes around this problem attacks the notion of "locality".
One possible way to resolve this is with complete observables. The idea is to consider gauge-invariant relations between gauge-dependent fields. The approach is motivated by the idea that all that matters in general relativity are relations between dynamical quantities. For a review of this, see Tamborino (arXiv:1109.0740) and for applications to general relativity (arXiv:gr-qc/0610060 and arXiv:gr-qc/0702093).
Another approach uses holonomies for observables (see J.N. Goldberg, J. Lewandowski and C. Stornaiolo,
  "Degeneracy in loop variables".
  Commun.Math.Phys. 148, no.377 (1992)
  doi:10.1007/BF02100867,
  Project Euclid's Eprint.). This uses a change of variables from the metric (and its conjugate momenta) to spin-connections (and its conjugate momenta) popularly called "Ashtekar variables", swapping the (local) Lorentz group $O(3,1)$ for the simpler $SL_{2}(\mathbb{C})$ group. These $SL_{2}(\mathbb{C})$-connections's holonomies are clearly nonlocal quantities, thus circumventing Torre's "no-go" theorem.
One last approach worth mentioning [arXiv:gr-qc/9708041 and arXiv:gr-qc/9906044] is to use the eigenvalues for the Dirac operator defined on a given spacetime. Although not all spacetimes admit a Dirac operator, the general idea may be generalized to the Klein-Gordon operator.
