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