the density field was resolved in a $1024^3$ resolution grid. Furthermore, the simulations were initialized at redshift $z_i=49$ using the Zel’dovich approximation [@Zeldovich:1969sb] and evolved through $z_f=0$.

Gravitationally bound dark matter halos were identified using the [ROCKSTAR]{} halo finder [@2013ApJ...762..109B]. Finally, so as to get an estimate of the variance, each model was simulated for 5 random realizations, corresponding to different random phases in the initial density field.

The $1024^3$ $(Mpc/h)^3$ volume simulations’ results in Group I become noisy at scales $r>100$ Mpc/h. To probe the BAO scales, where LPT has been previously found to perform very well for GR [@Vlah:2014nta], we also test our results against the largest volume $f(R)$ simulations performed to date for the modified gravity lightcone simulation project [@Arnold:2018nmv]. In these simulations, which we will call Group II from now on, the box side is $L_{box}=1536$ Mpc/h with $2048^3$ equal mass particles used, for GR and the $\abs{\bar{f}_{R_0}}$=$10^{-5}$ model. The parameters describing the background $\Lambda$CDM cosmology are the best-fit ones given by the Planck collaboration and have the following values: $\Omega_b=0.0486$, $\Omega_m=0.3089$, $\Omega_L=0.6911$, $h=0.6774$, $n_s=0.9667$ and $\sigma_8=0.8159$. The simulations were performed using the MG code [MG-GADGET]{} [@doi:10.1093/mnras/stt1575], which is a MG extension to the code [P-GADGET3]{}, 