What is the physical interpretation of harmonic coordinates? 
When I see harmonic coordinates used somewhere, what should my association be?
Is there some general use or need to consider the harmonic coordinate condition?

I don't really see what's behind all the things said in the wikipedia article and I've seen these coordinates used in some papers but don't really understand their significance, other than that one can consequently algebraically use the relation $\Delta x^\mu=0$ in computations. I don't see what that means though, do the full set of these coordinates have some general geometrically visible properties, maybe characterized with the associated base of distributions?
Moreover, do these coordinates have some special role in geometric quantization?
 A: The well known property of the harmonic coordinates is that the covariant
divergence of a vector field and the d'Alambertian of a scalar field take a particularly simple
form:
$$
D_{\mu}A^{\mu} \rightarrow g^{\mu\nu}\partial_{\mu}A_{\nu},\\
g^{\mu\nu}D_{\nu}D_{\mu}\phi \rightarrow g^{\mu\nu}\partial_{\mu}
\partial_{\nu}\phi.
$$
The harmonic condition
$$
\partial_{\mu}\left(  \sqrt{-g}g^{\mu\nu}\right)  =0\qquad\qquad\left(  1\right)
$$
is widely used to construct the so called de Donder gauge for the quantization of
a weak gravitational field. If one uses the deviation $\psi^{\mu\nu}$ of the
contravariant metric density from the flat one $\eta^{\mu\nu}=diag\left(
1,-1,-1,-1\right)  $ as a field variable:
$$
\sqrt{-g}g^{\mu\nu}=\eta^{\mu\nu}+\psi^{\mu\nu},
$$
then the gauge condition $\partial_{\mu}\psi^{\mu\nu}=0$ looks very similar to
the well known Lorentz condition (Feynman gauge) $\partial_{\mu}A^{\mu}=0$. If
one would like to use the deviation of the covariant metric tensor as a field
variable
$$
g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},\qquad\left(  2\right)
$$
then the weak field expansion of the condition (1) take the form:
$$
\partial_{\mu}\left(  h^{\mu\nu}-\frac{1}{2}\eta^{\mu\nu}h_{\alpha}^{\alpha
}\right)  =0.
$$
The weak-field expansion of the Einstein-Hilbert action with respect to
$h_{\mu\nu}$-field (2) has the form:
$$
S  =\frac{1}{16\pi G_{N}}\int\mathrm{d}^{4}x\,\sqrt{-g}\,R\\
 =\frac{1}{2\kappa^{2}}\int\mathrm{d}^{4}x\,\left[  \,\partial_{\alpha}
h_{\mu\nu}\partial^{\alpha}h^{\mu\nu}-\,\partial_{\alpha}h\,\partial^{\alpha
}h-2\,\partial_{\mu}h^{\mu\nu}\left(  \partial_{\alpha}h_{\nu}^{\alpha
}-\partial_{\nu}h\right)  +O\!\left(  h^{3}\right)  \right]  ,
$$
where $\kappa=\sqrt{32\pi G_{N}}$. The gauge can be fixed by adding the term:
$$
\frac{1}{\kappa^{2}}\int\mathrm{d}^{4}x\left(  \partial_{\alpha}h_{\mu
}^{\alpha}-\frac{1}{2}\partial_{\mu}h\right)  \left(  \partial_{\beta}
h^{\beta\mu}-\frac{1}{2}\partial^{\mu}h\right)  ,
$$
thus the action takes a particular simple form:
$$
S=\frac{1}{2\kappa^{2}}\int\mathrm{d}^{4}x\,\left[  \,\partial_{\alpha}
h_{\mu\nu}\partial^{\alpha}h^{\mu\nu}-\,\frac{1}{2}\partial_{\alpha
}h\,\partial^{\alpha}h+O\!\left(  h^{3}\right)  \right]  .
$$
Therefore, in the de Donder gauge the graviton propagator has a very simple form:
$$
D_{\mu\nu,\alpha\beta}=\left\langle 0\left\vert T\,h_{\mu\nu}\left(  x\right)
h_{\alpha\beta}\left(  y\right)  \right\vert 0\right\rangle =i\kappa^{2}
\int\frac{\mathrm{d}^{4}p}{\left(  2\pi\right)  ^{4}}\frac{e^{-ip\cdot\left(
x-y\right)  }}{p^{2}+i0}\times\frac{1}{2}\left(  \eta_{\mu\alpha}\eta_{\nu\beta
}+\eta_{\mu\beta}\eta_{\nu\alpha}-\eta_{\mu\nu}\eta_{\alpha\beta}\right)  .
$$
The de Donder gauge is so to speak the GR analogue of the Feynman gauge for QCD or QED.
Using the gauge condition (1) and vertices extracted from the weak-field
expansion of the Einstein-Hilbert action and utilizing the QFT perturbation theory
with respect to $h_{\mu\nu}$, one can find, for example, the gravitational
field of a static spinless source. The result will be no more than the $r_{g}/r$-expansion of the Schwarzschild metric in the harmonic coordinates (see, e.g., S. Weinberg, Gravitation and Cosmology, eq. (8.2.15)):
$$
ds^{2}=\frac{1-r_{g}/\left(  2r\right)  }{1+r_{g}/\left(  2r\right)  }
\,dt^{2}-\frac{1+r_{g}/\left(  2r\right)  }{1-r_{g}/\left(  2r\right)
}\,\frac{r_{g}^{2}}{4r^{4}}\left(  \mathbf{r}\cdot d\mathbf{r}\right)
^{2}-\left(  1+\frac{r_{g}}{2r}\right)  ^{2}d\mathbf{r}^{2}.
$$
A: One can use harmonic coordinates to define a Lorentz ether interpretation for the Einstein equations. In this interpretation, one of the harmonic coordinates should be time-like and used as absolute time, the other as the Euclidean coordinates of the absolute space background. Then, the harmonic conditions itself can be interpreted as the continuity and Euler equations of the Lorentz ether. 
$g^{00}\sqrt{-g} = \rho$, $g^{i0}\sqrt{-g} = \rho v^i$, $g^{ij}\sqrt{-g} = \rho v^i v^j - \sigma^{ij}$ are, then, the formulas which connect the basic properties of the Lorentz ether (density, velocity, stress tensor) with the gravitational field. 
All other fields have to define other properties of the ether. 
See http://ilja-schmelzer.de/ether
