Rigorous approach to Newtonian limit could be achieved within frame theory of J. Ehlers, which contains both Newtonian gravity and general relativity. Frame theory has a parameter, $\lambda$ which can be though of as $1/c^2$ (where $c$ is the speed of light). For $\lambda=0$ the frame theory coincides with Newton–Cartan theory, which in turn could be thought of as geometrized version of Newtonian gravity. For $\lambda>0$ the frame theory is just general relativity rewritten in a language modeled after Newton–Cartan theory.
Frame theory deals with the following collection of fields on a 4-manifold $M$, spacetime:
$t_{μν}$, a nowhere vanishing symmetric, 2-covariant tensor field, the temporal metric;
$s^{μν}$, a nowhere vanishing, symmetric, 2-contravariant tensor field, the (inverse) spatial metric;
$\Gamma^\alpha_{μν}$, a symmetric, linear connection, the gravitational field;
$T^{μν}$, a symmetric, 2-contravariant tensor field, the mass–momentum–stress, or matter tensor.
Newtonian limit limit in frame theory is then defined for a family of solutions of Einstein equations depending on $\lambda$ in such a way that the fields have a limit at $\lambda\to0$, consistent with equations of Newton–Cartan theory.
For example, for Minkowski spacetime in cartesian coordinates we have
$$t_{μν}=−λ η_{μν}=\mathop{\mathrm{diag}}(1,-λ,-λ,-λ), \qquad s^{μν}=η^{μν}=\mathop{\mathrm{diag}}(−λ,1,1,1).$$
For $\lambda=0$ these fields have a meaningful limit, describing empty space of Newton–Cartan theory.
Is there a mathematically rigorous justification of this such as schematically “given a small $ϵ$, starting from some initial data there exists a time depending on $ϵ$ until which an appropriate difference of the solution of the 0th order expansion and the actual evolution can be bounded by $ϵ$”?
Not exactly. One should keep in mind that initial data of general relativity must satisfy Einstein constraint equations and thus initial data for some $ϵ>0$ are generally would be incompatible with initial data of Newtonian case for $ϵ=0$. So for the Cauchy problem the algorithm would be:
Find the one-parametric family of initial data on spacelike hypersurface that solve Einstein constraint equations for $ϵ>0$ and converge in the limit $ϵ→0$ to initial data of Newtonian problem;
Once appropriate initial data are chosen, local existence theorems ensure that there exist one-parameter family of solutions for all $0<ϵ < ϵ_0$ in a spacetime region;
The nontrivial part is to demonstrate that there is a finite time interval for which solutions exist for all values of $ϵ \in (0,ϵ_0)$;
And that these solutions converge to a solution of Newtonian theory as $ϵ\to 0$, this convergence is uniform in $ϵ$ and once can show the existence of an order $ϵ$ error estimate for the difference between Newtonian and relativistic solutions.
If this could be achieved we could say that there exists a rigorous Newtonian limit.
There are indeed several results establishing rigorous Newtonian limits (in the sense outlined above) including various cosmologies, most notably in the works of Todd Oliynyk.
The application in my mind is in cosmology. How reliable is the Newtonian approximation in cosmology?
One should remember that while techniques outlined above do prove that there are GR solutions that are well described by Newtonian gravity this alone does not guarantee that the real world would satisfy the prerequisites under which these results were derived. As an example of issues that could arise in such situations one can check the debate on cosmological backreaction (a question not unrelated to the validity of Newtonian approximation): original papers (1, 2, 3, 4); criticism, response. Though mathematically results seem to be solid its applicability to the real world is questioned.
References
For general introduction to frame theory I would recommend a “Golden Oldies” republication (and translation from German) of 1981 paper by Ehlers:
Newton–Cartan cosmology:
Ehlers, J., & Buchert, T. (1997). Newtonian cosmology in Lagrangian formulation: foundations and perturbation theory. General Relativity and Gravitation, 29(6), 733-764, doi:10.1023/A:1018885922682, arXiv:astro-ph/9609036.
Ehlers, J. (1999). Some developments in Newtonian cosmology. In: Harvey A. (eds) “On Einstein’s Path” (pp. 189-202). Springer, New York, NY.
On Newtonian limits:
Cosmological Newtonian limit:
Oliynyk, T. A. (2014). Cosmological Newtonian limit. Physical Review D, 89(12), 124002,
doi:10.1103/PhysRevD.89.124002, arXiv:1307.6281.
Oliynyk, T. A. (2015). The Newtonian limit on cosmological scales. Communications in Mathematical Physics, 339(2), 455-512, doi:10.1007/s00220-015-2418-5, arXiv:1406.6104.
Liu, C., & Oliynyk, T. A. (2018). Cosmological Newtonian limits on large spacetime scales. Communications in Mathematical Physics, 364(3), 1195-1304, doi:10.1007/s00220-018-3214-9,
arXiv:1711.10896.
