$f(R)$ model [@Hu:2007nk], which we study in this paper and the functional form of which is given by: $$\label{fRHu}
f(R)=-m^2\frac{c_1\left(R/m^2\right)^n}{c_2\left(R/m^2\right)^n+1}.$$ In equation (\[fRHu\]), $\Omega_{m0}$ denotes the matter fractional energy density and $H_0$ the Hubble Constant, both evaluated today, $m=H_0\sqrt{\Omega_{m0}}$, which has dimensions of mass and $n$, $c_1$ and $c_2$ are free parameters.

In order to match the $\Lambda$CDM expansion history and for sufficiently small values of $|f_{R_0}|$, the background value of the Ricci scalar, $\bar{R}$, becomes equal to: $$\begin{aligned}
\bar{R}= 3 \Omega_m H_0^2 \left(1+ 4 \frac{\Omega_{\Lambda0}}{\Omega_{m0}} \right).\end{aligned}$$ where $\Omega_{\Lambda0}$ is the dark energy fractional density evaluated at the present time.

The derivative $f_{R}=\frac{df(R)}{dR}$ becomes a functional of the cosmological parameters when evaluated today, through the relationship $$\bar{f}_{R_0} =-n\frac{c_1}{c_2^2}\left(\frac{\Omega_{m0}}{3(\Omega_{m0}+\Omega_{\Lambda0})}\right)^{n+1}.$$ The above mapping allows us to reduce the number of free parameters and the Hu-Sawicki model is commonly parametrized by quoting the values chosen for $n$ and $|f_{R_0}|$, with the latter being the background value of the fifth potential evaluated today.

The reason the Hu-Sawicki model is so popular is that it can be cast into the form of a scalar-tensor theory that realizes the chameleon screening mechanism [@PhysRevD.69.044026; @PhysRevLett.93.171104], through a conformal transformation [@Brax:2008hh], with the quantity $f_{R}$ identified as the 