-1.08 -1.06 -1.04 -1.02 -1.00 -0.98 -0.96 -0.94 -0.92 0 1 2 3 4 w DE z n=2  =1 n=2  =3 n=2  =10 CDM -1.08 -1.06 -1.04 -1.02 -1.00 -0.98 -0.96 -0.94 -0.92 0 1 2 3 4 w DE z n=2  =0.95 n=3  =0.73 n=4  =0.61 CDM FIG. 1: Evolution of the equation-of-state parameter of effective dark energy. III. PERTURBATIONS We now turn to the evolution of density fluctuations. In f ( R ) gravity, the evolution equation of density fluctuations, δ , deeply in the sub-horizon regime is given by[7] In the previous paper[6] we obtained an analytic solution for the high-curvature regime when the scale factor evolves as a ( t ) ∝ t 2 / 3 and F takes an asymptotic form F ≃ 1 − 2 nλ ( R/R$_{s}$ ) − 2 n − $^{1}$. In the present paper, we numerically integrate (7) up to z = 0 without using the approximation | F − 1 | ≪ 1. The wavenumber of our particular interest is the scale corresponding to σ$_{8}$ normalization, for which we find k$_{eff}$ ( r = 8 h − $^{1}$Mpc) = 0 . 174 h Mpc − $^{1}$. Since the standard ΛCDM model normalized by CMB data explains galaxy clustering at small scales well, δ$_{fRG}$ should not be too much larger than δ$_{ΛCDM}$ on these scales. We may typically require ( δ$_{fRG}$/δ$_{ΛCDM}$ ) 2 ≲ 1 . 1. Although we neglect non-linear effects here, the difference between linear calculation and non-linear N-body simulation remained smaller than 5% at wavenumber 0 . 174 h Mpc − $^{1}$[8]. The left panel of Fig. 2 represents ( δ$_{fRG}$/δ$_{ΛCDM}$ ) 2 as a function of λ for n = 2 together with two fitting functions. The solid line is from the analytic formula obtained in Ref. [6], and the broken line is numerical fitting using an exponential function 1 + b$_{n}$e − q n $^{λ}$. From these analysis, we can constrain the parameter space as the right panel of Fig. 2. The region which satisfy ( δ$_{fRG}$/δ$_{ΛCDM}$ ) 2 < 1 . 1 lies above the solid line. The region below the dotted line 4 ¨ δ + 2 H ˙ δ − 4 πG$_{eff}$ρδ = 0 , G$_{eff}$ = G F 1 + 4 k 2 a 2 F ′ F 1 + 3 k 2 a 2 F ′ F . (7)