11 FIG. 3: (Colour Online) The gravitational potential on the z plane where z = 32 h − 1 Mpc. Dark regions are where the potential is more negative (deeper) while light regions are where it is less negative (shallower). The 4 columns are for the 4 models we consider and the 3 rows are for three output times, respectively a = 0 . 2 , 0 . 5 and 1 (corresponding to redshifts 4, 1 and 0) where a is the cosmic scale factor. should have remained in the halo would be incorrectly removed by MHF. In general, similar things should be taken care of in other theories involving modifications to gravity in the non-relativistic limit, such as MOND and f ( R ) gravity. An exact calculation of the escape velocity in the cou- pled scalar field model is obviously difficult due to the complicated behaviour of the scalar field, and thus here we introduce an approximated algorithm, which is based on the MHF default method [42], to estimate it. MHF works out v$_{esc}$ using the Newtonian result v 2 esc = 2 | Φ | in which Φ is the gravitational potential. Under the assumption of spherical symmetry, the Poisson equation ∇ $^{2}$Φ = 4 πGρ$_{m}$ could be integrated once to give suming spherical symmetry of the halo, calculate the es- cape velocity v$_{esc}$ at the position of each particle, (ii) if the velocity of the particle exceeds v$_{esc}$ then it does not belong to the virialized halo and is removed. (i) and (ii) are then iterated until all unbound particles are removed from the halo or the number of particles in the halo falls below a pre-defined threshold, which is 20 in our simula- tions. Note that the removal of unbound particles is not used in some halo finders using the spherical overdensity (SO) algorithm, which includes the particles in the halo as long as they are within the radius of a virial density contrast. Another advantage of MHF is that it does not require a pre-defined linking length in finding halos, such as the friend-of-friend procedure. When it comes to our coupled scalar field model, the MHF algorithm also needs to be modified. The reason is that, as we mentioned above, the scalar field ϕ behaves as an extra ”potential” (which produces the fifth force) and so the CDM particles experience a deeper total ”gravita- tional” potential than what they do in the ΛCDM model. Consequently, the escape velocity for CDM particles in- creases compared with the Newtonian prediction. This is indeed important to bear in mind because, as we have seen in Table I, the CDM particles are typically much faster than what they are in the ΛCDM simulation, and so if we underestimate v$_{esc}$ then some particles which which is just the Newtonian force law. This equation can be integrated once again to obtain where Φ$_{0}$ is an integration constant and can be fixed [42] d Φ dr = GM ( < r ) r 2 (39) Φ( r ) = G ∫ r 0 M ( < r ′ ) r ′ 2 dr ′ + Φ$_{0}$ (40)