0 ≤ z ≤ 1 to which all BAO and most of SN data refer is very close to − 1. Moreover, in this range (but not for larger values of z ), the behaviour of w$_{DE}$ for minimal allowed values of λ (i.e. for largest possible deviations from the ΛCDM background model) is well fitted by the CPL fit[25] w$_{DE}$ ( z ) = w$_{0}$ + w$_{a}$z/ (1 + z ) with ( n, λ$_{min}$, w$_{0}$, w$_{a}$ ) = (2 , 0 . 95 , − 0 . 92 , − 0 . 23), (3 , 0 . 73 , − 0 . 94 , − 0 . 22) and (4 , 0 . 61 , − 0 . 96 , − 0 . 21), respectively. | 1 + w$_{0}$ | and | w$_{a}$ | decrease slowly for larger values of n . These values of w$_{0}$ and w$_{a}$ lie very close to the center of the 68% and 95% CL ellipses for all combined data in Fig. 13 of Ref. [2]. As explained above, this phantom crossing behaviour is not peculiar to the specific choice of the function (3) but a generic one in models which satisfy the stability condition F ′′ > 0. Indeed, a similar behaviour has been observed in other f ( R ) DE models, too [18, 26]. We also note that different definitions of ρ$_{DE}$ , P$_{DE}$ , and w$_{DE}$ have been used in literature [27] which lead to different behaviour of w$_{DE}$ . Although the behaviour of dark energy is quite different depending on model parameters, the total expansion factor a$_{0}$/a$_{i}$ from the epoch Ω$_{i}$ = 0 . 998 to the present varies only between a$_{0}$/a$_{i}$ = 10 . 8 and 11, the latter corresponding to the value in the ΛCDM model. We have also calculated the quantity B ( z ) = ( f ′′ /f ′ )( dR/d ln H ) introduced in Ref. [28] at present time. We have found B (0) = 0 . 21, 6 . 1 × 10 − $^{5}$, and 0 . 17, for ( n, λ ) = (2 , 0 . 95), (2 , 8), and (4 , 0 . 61), respectively. We now turn to evolution of density fluctuations. In f ( R ) gravity, the evolution equation of density fluctuations, δ , deeply in the sub-horizon regime is given by [29, 30] III. DENSITY FLUCTUATIONS ¨ δ + 2 H ˙ δ − 4 πG$_{eff}$ρδ = 0 , (13) G$_{eff}$ = G F 1 + 4 k 2 a 2 F ′ F 1 + 3 k 2 a 2 F ′ F , F ( R ) ≡ f ′ ( R ) . (14) This equation reduces to the correct evolution equation for all wavenumbers for the CDM model in the Einstein gravity where F = 1. In the previous paper[15] we obtained an analytic solution in the high-curvature regime when the scale factor evolves as a ( t ) ∝ t 2 / 3 and F takes the asymptotic form F ≃ 1 − 2 nλ ( R R$_{s}$ ) − 2 $^{n}$− 1 ≡ 1 − ( R R$_{c}$ ) − $^{N}$− 1 , (15) 7 where