content_image/400102.jpg
<img src='content_image/400102.jpg'>

the same quantities at a = 1 . 0 (late times). We can see the trend of increasing deviation from ΛCDM for all 4 models, because fifth force is essentially unsuppressed at the late epoch [cf. Fig. 2]. For example, the deviation of the model γ = 1 . 0 , µ = 10 − 5 is significantly larger than that of the model γ = 0 . 5 , µ = 10 − 5 as na¨ıvely expected, thanks to the lack of suppression of fifth force in both models (the γ = 0 . 5 model obviously has a smaller satu- rated fifth force). 

For comparison we also plot the Δ P/P that is pre- dicted by the linear perturbation theory for the 4 mod- els under consideration (the dashed curves). As can be seen there, at large scales (small k ) where linear pertur- bation is considered as a good approximation, the lin- ear and nonlinear results agree pretty well (especially for the a = 0 . 5 case). The largest scale we can probe is limited by the size of our simulation boxes (64 Mpc/h ) and as a result we cannot make plot beyond the point k ∼ 0 . 1 h/Mpc , where nonlinearity is expected to first become significant. However, in the case of a = 1 . 0, we can see the clear trend of the linear and nonlinear results merging towards k < 0 . 1 h/Mpc at vanishing Δ P/P . Similar results can be found in Fig. 2 of [26] for f(R) gravity. 

Because we have two species of matter particles, one 

uncoupled to the scalar field, we are also interested in their respect power spectrum and the bias between them. These are displayed in Fig. 6. The results could be under- stood easily: because a CDM particle always feel stronger total force than a baryon at the same position, so the clustering of the former is identically stronger than the latter as well. This could result in a significant bias be- tween these two species at the present time, especially for the models with µ = 10 − $^{5}$, where the fifth force is less suppressed. 

We identify halos in our N -body simulations using MHF (MLAPM Halo Finder) [41], which is the default halo finder for MLAPM. MHF optimally utilizes the re- finement structure of the simulation grids to pin down the regions where potential halos reside and organize the refinement hierarchy into a tree structure. Because MLAPM refines grids according to the particle density on them, so the boundaries of the refinements are sim- ply isodensity contours. MHF collect the particles within these isodensity contours (as well as some particles out- side). It then performs the following operations: (i) as- 

## D. Mass Functions 

