In most treatments of General Relativity, when the the Einstein-Hilbert action over some manifold $mathcal{M}$ (plus Gibbons-Hawking-York term if $mathcal{M}$ has a boundary), given by

$$S=frac{1}{2kappa}int_{mathcal{M}}star,mathcal{R}+frac{1}{kappa}int_{partialmathcal{M}}star,K$$

(with $8pi G=kappa$) is varied, it is done so by an implicit choice of chart. Namely, one writes the action (locally) as

$$S=frac{1}{2kappa}intmathrm{d}^dx,sqrt{det(g)},mathcal{R}+frac{1}{kappa}intmathrm{d}^{d-1}y,sqrt{det(h)},K$$

From this, one uses the coordinate equations for $mathcal{R}$ and vary it with the components of the metric tensor in this particular chart to derive the Einstein field equations, which is a very long and tedious process, prone to mistakes from misplaced indices.

However, it is generally possible to derive equations of motion for field theories without referencing a coordinate chart. In fact, such methods are typically very generalizable and very powerful for calculating conserved quantities. Consider, for instance, (Euclidean) Maxwell theory of a $U(1)$ gauge field $A$. The action is simply

$$S=frac{1}{2e^2}int_{mathcal{M}}Fwedgestar,F+int_{mathcal{M}} Awedgestar,j,$$

where $F=mathrm{d}A$ and $j$ is some matter current. Under a general transformation $Ato A+delta A$, we have

$$delta S=frac{1}{e^2}int_{mathcal{M}}mathrm{d}delta Awedgestar,F+int_{mathcal{M}}delta Awedgestar,j\=int_{mathcal{M}}delta Awedgeleft(frac{1}{e^2}mathrm{d}star F+star,jright)+int_{partialmathcal{M}}delta Awedgestar,F,$$

from which the equations of motion are immediately $mathrm{d}F=0$

and $mathrm{d}star F=-e^2star,j$, assuming $delta A$ has no support on $partialmathcal{M}$.

Does there exist a similar completely covariant and coordinate-independent derivation of the Einstein field equations from the Einstein-Hilbert action?