Is Hilbert-Einstein action just the leading order of some kind of series? Introducing the action for the gravitational field my GR professor stated that, in principle, one could write it as
$$S = k\int d^4x\sqrt{g}(\sum_n\sum_m a_{nm} R_n^m - 2\Lambda), \space \space \space \space m,n\in \mathbb{N}\backslash \{ 0\}$$
where $R_n^m$ is the $m$-th power of the $n$-th independent scalar achievable from the Riemann tensor in four dimensions multiplied by a constant coefficient $a_{nm}$, $g$ is the determinant of the metric tensor, $\Lambda$ is the cosmological constant, and $k$ is an arbitrary constant; his claim was based on the need of $S$ to be a scalar, and the fact that $S$ should contain at least second derivatives of the metric tensor (not considering $\Lambda$). Assuming $R_1$ to be the Ricci scalar, the $m=n=1$ term of $S$ is the Hilbert-Einstein action $S_{HE}$.
He then observed that every term different from $S_{HE}$ depends, dimensionally, from $G^k$, $k\in \mathbb{N}$, assuming that $c=\hbar=1$ and $G$ is the only other fundamental constant; from this, he proceeded in eliminating this terms saying that at our scales we are not sensible to them, i.e, every gravitational system we are able to study experimentally so far obeys EFE as derived from $S_{HE}$.
Is this the correct way to pose the problem, i.e, is this equivalent to postulate the validity of EFE? And if so, does this implies that when studying a generic system using EFE we are just neglecting higher order terms?
 A: Whenever you talk about "validity" of a theory/equation in physics, you must also specify the length scales being probed. This is important because physics at different scales is different, with different degrees of freedom and different characteristics. This is why we don't need to worry about Heisenberg's uncertainty principle and the insanely complicated quantum dynamics of air molecules when designing an aircraft; Newton's laws and fluid dynamics work pretty well. This means that a theory $(T_A)$ that is used to explain phenomena at one particular length scale $(L_A)$ should not be used to explain phenomena at a different length scale $(L_B)$. We say that the theory $T_A$ is effective for length scales of the order $L_A$ while another theory $T_B$ is effective for length scales of the order $L_B$. In this sense, every physical theory is an effective field theory.
General relativity (GR) is an effective field theory too. The typical way to solve problems in physics is to do perturbation theory. If we take the Einstein-Hilbert action
$$ S_{EH} = \int d^4 x \ M_P^2 \sqrt{-g} R  \qquad \qquad M_P=\text{Planck mass}$$
and expand it around a background, say flat space ($\bar{g}_{\mu \nu} = \eta_{\mu \nu}$), using the perturbative expansion $g_{\mu \nu} = \bar{g}_{\mu \nu} + h_{\mu \nu}$, we get (schematically)
$$ S_{EH} = \int d^4 x \  M_P^2 \left(h \partial^2 h +h (h \partial^2 h) + h^2(h \partial^2 h) + \cdots  \right)  $$
where we have ignored total derivative terms. Indices are suppressed for clarity. The kinetic term should not have coupling, so canonically normalizing the graviton field $h_{\mu \nu} \to \frac{h_{\mu \nu}}{M_P}$ gives an infinite series of graviton-graviton vertices suppressed by appropriate powers of $M_P$:
$$ S_{EH} = \int d^4 x \   \left(h \partial^2 h + \frac{h}{M_P} (h \partial^2 h) + \frac{h^2}{M_P^2} (h \partial^2 h) + \cdots  \right)  $$
Note the crucial negative power of coupling $M_P$ in interactions. Fourier transforming $\partial^2 \to -k^2$, we see that the perturbative expansion works well for energy scales $k \ll M_P$. As we probe extremely short length scales around Planck $\sim 10^{-35} m$, momenta scale like $k \sim M_P$, and all interaction terms in the perturbative expansion become equally important. The theory becomes strongly coupled and leads to a breakdown of perturbation theory. GR becomes non-perturbative for $k \gtrsim M_P$. If we use just perturbation theory, we don't know what happens to GR for Planck scale energies: GR is non-renormalizable. Actually, pure GR is renormalizable at one loop (in four dimensions, because of Gauss-Bonnet identity) as was shown by 't Hooft and Veltman. They showed, however, that gravity coupled to scalar matter is non-renormalizable at one loop with the counterterm Lagrangian
$$ \Delta \mathcal{L} \sim \frac{1}{\epsilon} \sqrt{-g} \ R^2$$
Goroff and Sagnotti finally showed that pure GR is non-renormalizable at two loops, where the counterterm Lagrangian goes like
$$  \Delta \mathcal{L} \sim \frac{1}{\epsilon} \sqrt{-g} \ R_{ab}^{\ \ c d} R_{cd}^{\ \ ef} R_{ef}^{\ \ ab}$$
The original Eisntein-Hilbert Lagrangian doesn't contain these terms: GR is non-renormalizable. But these higher derivative curvatures force us to include them in the original Lagrangian if we are to make perturbative (also ambitiously unitary) sense of quantum field theory for gravity at Planck length scales. This is also step $1$ of constructing an effective field theory (EFT): given a Lagrangian, add all possible terms consistent with symmetries of the theory and arrange them in an energy expansion, that is, a sum of higher and higher dimensional operators, suppressed by higher powers of some cutoff scale $\Lambda$. $\Lambda$ is where the EFT breaks down and you need a UV completion (string theory?) to go beyond.
Taking these ideas into consideration, we construct, by hand, an effective field theory of gravity with the action
$$ S = \int d^4 x \sqrt{-g} \ M_P^2 R + c_1 R^2 + c_2 R_{\mu \nu} R^{\mu \nu} + c_3 R_{\mu \nu \rho \sigma} R^{\mu \nu \rho \sigma} + \text{higher order terms} $$
Note that $c_1, c_2, c_3$ are dimensionless and go like $\sim \frac{M_P^2}{M_s^2}$. $M_s$ is the energy scale where these higher derivative corrections become important. Doing again a perturbative expansion and canonically normalizing, the Einstein-Hilbert term in (Fourier space) goes like $k^2$ suppressed by higher and higher powers of $M_P$. The higher curvature terms go like $k^4$ suppressed by higher and higher powers of $M_P$ and $M_s$. For low enough momenta $k \ll M_s$, these higher curvature terms are severely suppressed and not important compared to the Einstein-Hilbert term. So although we expect these higher curvature terms to be there (also predicted by string theory), they don't affect our daily life on earth. For length scales on earth, the Newtonian potential is enough to explain physics. But for black holes and the singularities at the beginning of the universe, it's been a puzzle for decades now. And perhaps many more to come.
Addendum: Perhaps, the one higher curvature correction that makes an interesting appearance which is useful not just for theoretical studies but also for cosmological observations is the $R^2$ term in the form of Starobinsky inflation.
