explained in Section \[sec:Models\], an additional degree of freedom is present, that directly couples to matter and causes particles to deviate from the nominal geodesics of GR. Consequently, equations (\[eq:geopoisson\]) are, in principle, modified for a scalar-tensor theory and so is the LPT framework presented above. In [@Valogiannis:2016ane], the LPT approach was expanded for chameleons and symmetrons, including the first order contribution to the Klein-Gordon equation and was shown to perform very well in the context of the COLA hybrid framework. The LPT approach for scalar-tensor theories up to third order was presented in [@Aviles:2017aor], the main results of which we summarize in this section.

In [@PhysRevD.92.063522], which developed an SPT framework for studying MG theories with a screening mechanism in the nonlinear regime, based on the closure theory approximation in [@Taruya:2007xy], the scalar field KG equation for a Brans-Dicke-like (BD) theory with interactions of a scalar field, $\phi$, was written as: $$\label{BransDicke}
(3 + 2 \omega_{BD})\frac{1}{a^2} k_x^2\phi(\bold{k}_x,t) = 8 \pi G \bar{\rho}_m \delta(\bold{k}_x,t) - \mathcal{I}(\phi),$$ where $\mathcal{I}(\phi)$ denotes the perturbative form of the field self-interaction term: $$\begin{aligned}
\label{interaction}
 \mathcal{I}(\phi) &=& M_1(\bold{k},t)\phi \nonumber
\\
&& + \frac{1}{2}\int \frac{d^3 k_1 d^3 k_2}{\left(2 \pi\right)^3} \delta_D(\bold{k}-\bold{k}_{12})M_2(\bold{k}_{1},\bold{k}_{2})\phi(\bold{k}_{1})\phi(\bold{k}_{2})\nonumber\\
&& + \frac{1}{6}\int \frac{d^3 k_1 d^3 k_2 