In Fig. 1 we have shown some snapshots of the distri- bution of baryonic and CDM particles, to give some idea about the hierachical structure formation and for com- parisons with other figures below. It shows clearly how some clustering objects develop, with filaments connect- ing them together. The baryons roughly follow the clus- tering of CDM particles, but in some low density regions they become slightly separated. To understand how the motion of the particles is al- tered by the coupling to the scalar field, in Fig. 2 we have shown the correlation between the magnitudes of the fifth force and gravity on the CDM particles (remember that baryons do not feel the fifth force). Ref. [23] has made a detailed qualitatively analysis about the general trend of this correction, and here we just give a brief description: From Eqs. (31, 32) we could see that, when the scalar field potential , i.e. , the last term of Eqs. (31, 32) could be neglected , then the scalar field ϕ is simply proportional to the gravitational potential Φ and as a result Eq. (30) tells us that the strength of the fifth force is just 2 γ 2 times that of gravity; in other words, the effect of the scalar field is a rescaling of the gravitational constant by 1 + 2 γ $^{2}$. This is because in this situation, the effective mass of the scalar field, which is given by m 2 eff = ∂ $^{2}$V$_{eff}$ /∂ϕ $^{2}$, where is the effective total potential, is light and the fifth force is long-range, like gravity. For comparison, in Fig. 2 we also plot this 2 γ 2 proportion between the two forces, as a straight line: lg( f ) = lg(0 . 8 g )+lg(2 γ $^{2}$), where f, g denote respectively the magnitudes of the fifth force and gravity, and the factor 0.8 in front of g comes from the fact that only 80% of the particles are CDM (and thus contribute to the fifth force). This scaling relation actually sets an upper limit on how strong the fifth force could be relative to gravity, should it not be suppressed by other effects. In contrast, when the value of ϕ is small, the last term of Eqs. (31, 32) is not negligible and the scalar field ac- quires a heavy mass, making it short-ranged. As a result, a particle outside a high density region might not feel the fifth force exerted by particles in that region, even it is quite close to the region. But because it can feel gravity from that region, so the total fifth force on the particle becomes less than the 2 γ 2 scaling. In general, the value of ϕ is determined by µ, γ, ρ$_{CDM}$ as well as its background value ¯ ϕ (which sets the bound- ary condition to solve the interior value). At early times ¯ ϕ is very close to 0 and ρ$_{CDM}$ is high everywhere, making ϕ small everywhere too, and suppressing the fifth force so that it is significantly below the 2 γ $^{2}$-scaling (first row of Fig. 2). At later times, ¯ ϕ increases and ρ$_{CDM}$ decreases, weakening the above effect so that the fifth force becomes ”saturated” ( i.e. , approaches the 2 γ 2 prediction) and the points in the figure hit the straight lines (last two rows). Because decreasing µ and increasing γ have the same ef- fects of making ϕ small, in the models with µ = 10 − 6 the fifth force saturates later than in the models with µ = 10 − $^{5}$. In addition, because high ρ$_{CDM}$ tends to de- The lower solid curves in each panel of Fig. 5 display crease ϕ and increase the scalar field mass, so in high density regions (where gravity is stronger) the fifth force also saturates later. The agreement between the numerical solution of the fifth force and the 2 γ $^{2}$-scaling relation in cases of weak- chameleon effect serves as an independent check of our numerical code. Fig. 3 plots the spatial configuration for the gravita- tional potential at the same position and output times as in Fig. 1. As expected, the potential is significantly deeper where there is significant clustering of matter [cf. Fig. 1]. We also show in Fig. 4 the spatial configuration for the scalar field ϕ at the same output position and times. At early times when ϕ is small and the scalar field mass is heavy, the fifth force is so short-ranged that ϕ only de- pends on the local density. This means that the spatial configuration of ϕ in this situation could well reflect the underlying dark matter distribution, a fact which could be seen clearly in the first row. As time passes by, the mass of the scalar field decreases on average and ¯ ϕ in- creases, the value of ϕ at one point is more and more influenced by the matter distribution in neighboring re- gions, and such an ”averaging” effect weakens the con- trast and makes the plots blurring (last two rows). Fur- thermore, obviously decreasing µ and increasing γ could increase the scalar field’s mass, shorten the range of the fifth force, make ϕ less dependent on its value in neigh- boring regions, and thus strengthen the contrast in the figures. To have a more quantitative description about how the matter clustering property is modified with respect to the ΛCDM model, we consider the matter power spectra P ( k ) in our simulation boxes. The nonlinear matter power spectrum in the present work is measured using POWMES [40], which is a public available code based on the Taylor expansion of trigono- metric functions and yields Fourier modes from a number of fast Fourier transforms controlled by the order of the expansion. We also average the results from the 5 runs for each model and calculate the variance. In Fig. 5 shown are the fractional differences between the P ( k ) for our 4 models and for ΛCDM, at two different output times. At early times ( a = 0 . 5, upper solid curves in each panel), the difference is generally small, but still the 2 models with µ = 10 − 5 show up to ∼ 25% deviation from ΛCDM prediction. This is because for larger µ the scalar field is lighter and the fifth force less suppressed, its influence in the structure formation therefore enhanced. Notice that on small scales the deviation from ΛCDM de- creases, which is a desirable property of chameleon mod- els which are designed to suppress the fifth force on small scale high density regions. C. Power Spectrum V$_{eff}$ ( ϕ ) = V ( ϕ ) + C ( ϕ ) ρ$_{CDM}$ (38) 9