d^3 k_3}{\left(2 \pi\right)^6} \delta_D(\bold{k}-\bold{k}_{123})\nonumber \\&& \   \hspace{1.5cm} \times M_3(\bold{k}_{1},\bold{k}_{2},\bold{k}_{3})\phi(\bold{k}_{1})\phi(\bold{k}_{2})\phi(\bold{k}_{3}),\end{aligned}$$ where we adopted the standard notation $\bold{k}_{ijk}=\bold{k}_{i}+\bold{k}_{j}+\bold{k}_{k}$ and $M_1(\bold{k},t)$, $M_2(\bold{k}_{1},\bold{k}_{2})$ and $M_3(\bold{k}_{1},\bold{k}_{2},\bold{k}_{3})$ are mass terms. The higher order piece in the Fourier space representation of the interaction term in equation (\[interaction\]), incorporates the screening effect, up to third order, that is responsible for recovering GR at small scales. The mapping between a given scalar-tensor theory and the BD form above can be easily performed through assigning appropriate values to the mass terms and the BD coupling $\omega_{BD}$ above, as we will later show for our two models of study.

The perturbed modified Einstein equations have the form [@Aviles:2017aor] $$\begin{aligned}\label{eq:modpoisson}
\nabla_{\bold{x}}\mathcal{\hat{T}}\bold{\Psi} &= -\frac{1}{a^2}\nabla^2_{\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) -\frac{1}{2 a^2} \nabla^2\phi - \frac{1}{2 a^2}\left(\nabla^2_{\bold{x}}\phi -\nabla^2\phi \right),
\end{aligned}$$ where we use (\[Lagpos\]) and introduce the time derivative operator $\mathcal{\hat{T}}=\frac{d^2}{dt^2} + 2H\frac{d}{dt}$, as in [@Matsubara:2007wj]. The last term in the second line of (\[eq:modpoisson\]), called $\it{frame}$-$\it{lagging}$ in [@Aviles:2017aor], is a geometrical term that occurs due to the fact that, in LPT, the KG equation should be expressed in Lagrangian $\left(\nabla^2\right)$, rather than Eulerian coordinates $\left(\nabla^2_{\bold{x}}\right)$. Taking this into account, equations (\[BransDicke\])-(\[eq:modpoisson\]) are combined to give $$\begin{aligned}\label{eq:finpsi}
& \left(J^{-1}\right)_{ij}\mathcal{\hat{T}}\Psi_{i,j}(\bold{k}) = -A(k)\delta(\bold{k}) + \frac{k^2}{a^2 