Assuming $t_{0}$ refers to an epoch early enough that the density perturbations around the background density $\bar{\rho}$ are negligible, meaning $\rho_{m}(\bold{q},t_0)=\bar{\rho}_m$, allows us to obtain the dark matter fractional overdensity, $\delta_m$, in the Lagrangian picture: $$\label{delJac}
1 + \delta_m(\bold{x},t) = \int d^3q \delta_{D}\left[\bold{x}-\bold{q}-\bold{\Psi}(\bold{q},t)\right] = \frac{1}{J(\bold{q},t)},$$ with $\delta_{D}$ being the Dirac delta function and $J(\bold{q},t)$ the determinant of the deformation matrix $$\label{Jacobian}
J_{ij} = \frac{\partial x^i}{\partial q^j} = \delta_{ij} + \frac{\partial \Psi^i}{\partial q^j}.$$ For an irrotational flow, which is a good approximation for cold dark matter and assuming that the gravitational evolution is governed by GR, perturbations around a flat FRW metric give the geodesic and Poisson equations, in the quasi-static approximation and for sub-horizon scales, as: $$\begin{aligned}\label{eq:geopoisson}
\ddot{\bold{x}} + 2H\dot{\bold{x}} &= -\frac{1}{a^2}\nabla_{\bold{x}}\psi(\bold{x},t), \\
\frac{1}{a^2}\nabla^2_{\bold{x}}\psi(\bold{x},t) &=  4 \pi G\bar{\rho}_m \delta(\bold{x},t).
\end{aligned}$$ We should point out that in (\[eq:geopoisson\]) $\psi(\bold{x},t)$ denotes the metric perturbation, which should not be confused with the Lagrangian displacement field $\bold{\Psi}(\bold{q},t)$. In the LPT picture, we perturbatively expand $\bold{\Psi}$ as $$\label{eq:psiexp}
\bold{\Psi}(\bold{q},t) = \sum_{n=1}^{\infty}\bold{\Psi}^{(n)}(\bold{q},t) = \bold{\Psi}^{(1)}(\bold{q},t)+\bold{\Psi}^{(2)}(\bold{q},t)+\bold{\Psi}^{(3)}(\bold{q},t)...$$ and equations (\[Lagpos\])-(\[eq:psiexp\]) form a closed system, that is recursively solved for the various orders of $\bold{\Psi}$. The first order solution is the so-called Zel’dovich approximation [@Zeldovich:1969sb].

In MG theories, as also 