\mathcal{\hat{T}}D^{(1)}(k,t)=A(k)D^{(1)}(k,t).$$ In the GR limit, $A(k)=A(k=0)=4 \pi G\bar{\rho}_m$, $D^{(1)}$ becomes scale independent and is nothing else than the first order growing mode for GR; the Zel’dovich approximation.

Moving on to the second order piece, we have: $$\label{eq:Zeldisp2}
\bold{k}\cdot\bold{\Psi}^{(2)} = \frac{i}{2} \int \frac{d^3 k_1 d^3 k_2}{\left(2 \pi\right)^3} \delta_D(\bold{k}-\bold{k}_{12}) D^{(2)}(\bold{k}_1,\bold{k}_2)\delta^{(1)}_1\delta^{(1)}_2,$$ where, for compactness, we adopted the notation $\delta^{(1)}_1=(\bold{k}_1,0)$ and the second order growth factor, $D^{(2)}(\bold{k}_1,\bold{k}_2)$, is given by [@Aviles:2017aor]: $$\begin{aligned}
\label{eq:secgrowth}
D^{(2)}(\bold{k}_1,\bold{k}_2) &=& D^{(2)}_a(\bold{k}_1,\bold{k}_2)-D^{(2)}_b(\bold{k}_1,\bold{k}_2)\frac{\bold{k}_1\cdot\bold{k}_2}{k^2_1k^2_1} \nonumber
\\ && +D^{(2)}_{FL}(\bold{k}_1,\bold{k}_2)-D^{(2)}_{\delta \mathcal{I}}(\bold{k}_1,\bold{k}_2).\end{aligned}$$ The four individual components are given by: $$\begin{aligned}\label{eq:D2sources}
\left(\mathcal{\hat{T}}-A(k)\right)D^{(2)}_a(\bold{k}_1,\bold{k}_2) &=  A(k)D^{(1)}(k_1)D^{(1)}(k_2), \\
\left(\mathcal{\hat{T}}-A(k)\right)D^{(2)}_b(\bold{k}_1,\bold{k}_2) &=  \left(A(k_1)+A(k_2)-A(k)\right)D^{(1)}(k_1)D^{(1)}(k_2), \\
\left(\mathcal{\hat{T}}-A(k)\right)D^{(2)}_{\delta \mathcal{I}}(\bold{k}_1,\bold{k}_2) &=  \left(\frac{2A_0}{3}\right)^2\frac{k^2}{a^2}\frac{M_{2}(\bold{k}_1,\bold{k}_2)D^{(1)}(k_1)D^{(1)}(k_2)}{6\Pi(k)\Pi(k_1)\Pi(k_2)}, \\
\left(\mathcal{\hat{T}}-A(k)\right)D^{(2)}_{FL}(\bold{k}_1,\bold{k}_2) &=  \left(\frac{M_1}{3\Pi(k)}\right)K^{(2)}_{FL}(\bold{k}_1,\bold{k}_2)D^{(1)}(k_1)D^{(1)}(k_2).
\end{aligned}$$ The two last terms represent the second order contributions to the growth factor, given by the screening and frame-lagging effects, correspondingly, while the expression for $K^{(2)}_{FL}$ is given in [@Aviles:2017aor]. Despite its lengthier expression, when taking the GR limit we get $D^{(2)}_{FL}=D^{(2)}_{\delta \mathcal{I}}=0$ and $D^{(2)}_{a}=D^{(2)}_{b}$, allowing $D^{(2)}$ to become scale-independent, reducing to the known GR result, which can be well approximated by $D^{(2)}(t)=-\frac{3}{7}\left(D^{(1)}(t)\right)^2$ for $\Lambda$CDM cosmologies [@Bouchet:1994xp].

Solving for the third order piece in (\[eq:finpsi\]) results in a lengthy differential equation for the third order MG growth factor, $D^{(3)}(\bold{k}_1,\bold{k}_2,\bold{k}_3)$, that also needs to be symmetrized. The result is 