Can Einstein-Hilbert action be derived from symmetry considerations? The action of a free relativistic classical field theory can be derived from Poincare invariance, locality, and retaining terms quadratic in fields. Is there a similar set of symmetry principles which enable one to write down the Einstein-Hilbert (EH) action $$S_{EH}=\frac{c^4}{16\pi G}\int d^4x \sqrt{-g}R\tag{1}$$ uniquely?
In chapter VII.1, of the QFT book by A. Zee, it is mentioned that

The Einstein-Hilbert action is uniquely determined if we require the action to be coordinate invariant and to involve two powers of spacetime derivative.

which is then justified a posteriori. But how does one arrive at the form (1) from scratch and in a neat step-by-step manner?
 A: Well, one can reason as follows:


*

*One wants a diffeomorphism-invariant action, which must be of the form $$ S=\int d^4x\sqrt{-g}L, $$ where $L$ is a scalar in terms of transformation properties.

*Gravity is assumed to be purely metric, so $g_{\mu\nu}$ is the only dynamical variable that can appear in the action.

*So we must look for scalars that can be constructed out of the metric. One is of course $g^{\mu\nu}g_{\mu\nu}=4$, but since we can rescale by constants, this is basically a constant. Normally you wouldn't care about constants, but constants here do contribute, because they get multiplied by the metric determinant.

*At an arbitrary point $p$, you can set up coordinates such that $g_{\mu\nu}(p)=\eta_{\mu\nu}$ and $\partial_\mu g_{\nu\rho}(p)=0$, so there are no scalars that can be constructed out of the metric and its first derivatives. One needs to look for higher derivatives.

*The curvature tensor $R^\rho_{\ \sigma\mu\nu}$ is a good candidate. There are ways to show that the Riemann-tensor is basically unique or "most general" such a tensor, so we don't really need to look futher.

*It should be noted that one reason we don't want the Lagrangian to contain second or higher derivatives is due to Ostrogradsky instabilities. Ostrogradsky instability, is however avoided, if the Lagrangian is linear in the second derivatives.

*It is easy to check that the Riemann curvature tensor is linear in the second derivatives of the metric components.

*One needs to form a scalar out of $R^\rho_{\ \sigma\mu\nu}$, however any scalar produced by squaring a curvature tensor will ruin the linearity in the second derivatives. This rules out $R^{\rho\sigma\mu\nu}R_{\rho\sigma\mu\nu},\ C^{\rho\sigma\mu\nu}C_{\rho\sigma\mu\nu},\ R_{\mu\nu}R^{\mu\nu}$ etc. However the curvature scalar $R=R^\rho_{\mu\rho\nu}g^{\mu\nu}$ is still linear in the second derivatives.

*This gives us $$ L=\alpha R+\beta $$ as a possible scalar Lagrangian, where $\alpha$ and $\beta$ are constants. Of course $\beta$ may be zero, $\alpha$ not so much. We factorize as $$ L=\alpha\left(R+\frac{\beta}{\alpha}\right), $$ with $\alpha=\frac{1}{2\kappa}$ and $\beta=-\frac{\Lambda}{\kappa}$ to get $$ L=\frac{1}{2\kappa}\left(R-2\Lambda\right), $$ which is the Lagrangian for the Einstein Field Equation (extended with the cosmological constant), provided that $\kappa=\frac{8\pi G}{c^4}$.
