What is the physical interpretation of harmonic coordinates?

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}. $$