In what sense exactly are the Einstein metrics distinguished? EDIT: In general relativity given a manifold $M$ one can consider a functional on (pseudo-) Riemannian metrics $g$ $$\int_M R\,\, dvol_g,$$
where $R$ is the scalar curvature and $vol_g$ is the (pseudo-) Riemannian measure. The extremal metrics for this functional (solutions of the Euler-Lagrange equation) satisfy the so called Einstein equation.
As far as I heard the distinguished property of the above functional is that $R$ involves second derivatives of the metric, but the Euler-Lagrange equation is a differential equation not of the fourth order in the metric, but only of the second order in the metric.

I am looking for a precise statement within what class of functionals on metrics the above property distinguishes the scalar curvature $R$. Say assume I have a generally covariant expression involving at most second order derivatives of the metric such that the Euler-Lagrange equation is also of second order in the metric. Does it imply that this expression is proportional to $R$?

 A: If I understood your question correctly, the answer indeed is due to Lovelock. I think it's important to state all the hypotheses clearly, because they are not always reported accurately.
Theorem. (Lovelock, 1971) Given a metric $g_{ab}$ and a covariantly constructed symmetric 2-tensor $T^{ab}(g,\partial g, \partial^2 g)$ such that $\nabla_a T^{ab} = 0$, then necessarily the tensor $T^{ab} dvol_g = \frac{\delta L}{\delta g_{ab}}$ is the Euler-Lagrange derivative of the Lagrangian density $L = \sum_k \alpha_k R_{[a_1 a_2}{}^{a_1 a_2} \cdots R_{a_{2k-1} a_{2k}]}{}^{a_{2k-1} a_{2k}} dvol_g$, where $R_{abcd}$ is the Riemann curvature tensor, $[{\cdots}]$ denotes full antisymmetrization, and $\alpha_k$ are arbitrary constants (their values determine $T^{ab}$).
The number of terms in $L$ is finite in any given dimensions, because antisymmetrizing over more indices than the dimension always gives zero. In 4 dimensions, the only possibilities are $L = (\alpha_0 + \alpha_1 R + \alpha_2 R_{abcd} R^{abcd}) dvol_g$.
If $T^{ab} dvol_g = \delta L/\delta g_{ab}$, for any covariantly constructed Lagrangian $L$, then $T^{ab}$ is obviously symmetric and $\nabla_a T^{ab}=0$. So Lovelock's theorem actually classifies all Lagrangians that depend depend covariantly on the metric and Riemann curvature (including its derivatives) that have second order Euler-Lagrange equations, up to null Lagrangians, those whose Euler-Lagrange equations are identically zero. In fact, in 4 dimensions, a particular choice of the constants with $\alpha_2\ne 0$ gives such a null Lagrangian, the so-called Gauss-Bonnet term. So, among 4-dimensional Lagrangians, by adding a null Lagrangian, one can always reduce the desired variational principle to the Einstein-Hilbert form, $L= (\alpha_0 + \alpha_1 R) dvol_g$.

Lovelock, D., The Einstein tensor and its generalizations, J. Math. Phys. 12, 498-501 (1971). ZBL0213.48801.

A: OK I think I understand the question now.
You want to look at Lovelock gravities.
These are the most general theories (not bothering with a list of qualifiers here, see wikipedia) that produce 2nd order EOMS.
For dimensions d=3 and d=4 the only Lovelock functional is indeed the Einstein-Hilbert one.
