3\Pi(k)}\delta\mathcal{I}(\bold{k}) \\
& + \frac{M_1(k)}{3\Pi(k)}\frac{1}{2 a^2}\left(\nabla^2_{\bold{x}}\phi -\nabla^2\phi \right)(\bold{k}),
\end{aligned}$$ where, following the definitions in [@PhysRevD.92.063522], we have $$\begin{aligned}\label{eq:sources}
A(k) &=  4 \pi G\bar{\rho}_m\left(1 + \frac{k^2}{a^2 3\Pi(k)}\right), \\
\Pi(k) &= \frac{1}{3a^2}\left[\left(3+2\omega_{BD}\right)k^2 + M_1a^2\right]
\end{aligned}$$ and all the quantities are Fourier transforms in the Lagrangian $\bold{q}$-space. The inverse Jacobean in (\[eq:finpsi\]) reflects the derivative transformation to the $\bold{q}$-space, where the Einstein notation is adopted. Furthermore, $\delta(\bold{k})$, $\delta\mathcal{I}(\bold{k})$ and $\left(\nabla^2_{\bold{x}}\phi -\nabla^2\phi \right)(\bold{k})$ are the Lagragian Fourier transformations of the Lagrangian-transformed overdensity (\[delJac\]), the higher order interaction Kernels in (\[interaction\]) and the frame-lagging Kernel, correspondingly. The expression for the latter is given in [@Aviles:2017aor]. Equation (\[eq:finpsi\]) forms a closed system with (\[delJac\]), (\[Jacobian\]), (\[BransDicke\]) and (\[interaction\]) that is solved, perturbatively, to obtain the MG solution up to various orders in $\bold{\Psi}$, as in (\[eq:psiexp\]).

Solving for the first order solution, one gets [@Valogiannis:2016ane]: $$\label{eq:Zeldisp}
\bold{k}\cdot\bold{\Psi}^{(1)} = i D^{(1)}(k,t)\delta^{(1)}(\bold{k},t=0),$$ which can be easily solved for the displacement field, as: $$\label{eq:Zeldispsol}
{\Psi}^{j}(k,t) = \frac{i k^j}{k^2} D^{(1)}(k,t)\delta^{(1)}(\bold{k},t=0).$$ We see that the r.h.s of (\[eq:Zeldispsol\]) can be conveniently decomposed into a product of the first order density mode at very early times, $\delta^{(1)}(\bold{k},t=0)$, early enough to be gaussian and a space-time dependent growth factor $D^{(1)}(k,t)$, given by: $$\label{growth1st}
