scalar field that is coupled to matter. Through the interplay between matter and the self-interaction potential, the scalar chameleon field becomes massive near high over-densities and the associated fifth forces get exponentially suppressed due to the Yukawa effect. A stronger screening effect is manifested in lower values of the parameter $|f_{R_0}|$, with the limit of $|f_{R_0}|\rightarrow0$ exactly recovering GR. Following the literature and also because of the available simulations (as we will discuss below in \[sim\]), we choose $n=1$ and consider three different $f(R)$ models with $\abs{\bar{f}_{R_0}}=\{10^{-6},10^{-5}, 10^{-4}\}$, which we shall refer to, from now on, as F6, F5 and F4.

When considering perturbations around a homogeneous and isotropic FRW metric in the conformal Newtonian gauge, the resulting system of the Poisson and KG equations becomes [@Hu:2007nk]: $$\begin{aligned}\label{eq:poisson}
\nabla^2 \Phi_N &= 4 \pi G  a^2 \delta \rho_m - \frac{1}{2} \nabla^2 f_R, \\
\nabla^2 f_R &= -\frac{a^2}{3} \delta R - \frac{8 \pi G a^2}{3} \delta \rho_m,
\end{aligned}$$ where $\Phi_N$ is the Newtonian potential and $\delta \rho_m$ the matter density perturbation. $\delta R$, the perturbation to the Ricci scalar, can be written as a function of the scalar field $f_R$, which will play a central role in our screening implementation in Section 