Linearizing Gravity to ${\cal O}(h^3)$

Most of the difficulties of these calculations comes from the fact that people write the Lagrangian directly in terms of the metric perturbation $h_{\mu\nu}$.

It is much simpler to write it in terms of the difference of connections tensor $F_{\mu\nu}{}^\beta$, i.e., $$(\nabla_\mu-\bar{\nabla}_\mu)A_\nu = \mathcal{F}_{\mu\nu}{}^\beta{}A_\beta,$$ where $\nabla_\mu$ is the connection compatible with the full metric ($g_{\mu\nu}$) and $\bar{\nabla}_\mu$ compatible with the background metric ($\bar{g}_{\mu\nu}$ in your case Minkowski). The curvature scalar is naturally quadratic in $\mathcal{F}_{\mu\nu}{}^\beta$, hence, the higher order terms comes from the expansion of $\sqrt{-g}$ and $\delta{}g^{\mu\nu}$ in terms of the metric difference $\xi_{\mu\nu} = g_{\mu\nu} - \bar{g}_{\mu\nu}$. Using this method one don't need to choose a gauge a priori.

In our paper http://arxiv.org/abs/1206.4374 we develop the Lagrangian and Hamiltonian up to second order using the methodology described above. Additionally, we write the action in term of the kinetic quantities defined in geodesic space-time foliation. Using the results of this paper it is easy to generalize (in the Minkowski background) the Lagrangian to higher orders.

I think the earliest papers where this was written down were by DeWitt. But for a reference easily available via the arXiv look at hep-th/9411092. Eq (2.17) has the expansion you want and eq. (2.18) even has it to fourth order in $h$.