{'fraction_non_alphanumeric': 0.04527904527904528, 'fraction_numerical': 0.018252018252018252, 'mean_word_length': 4.374646293152235, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Can f(R) Modified Gravity Theories Mimic a LCDM Cosmology?

We consider f(R) modified gravity theories in the metric variation formalism and attempt to reconstruct the function f(R) by demanding a background LCDM cosmology. In particular we impose the following requirements: a. A background cosmic history H(z) provided by the usual flat LCDM parametrization though the radiation (w_eff=1/3), matter (w_eff=0) and deSitter (w_eff=-1) eras. b. Matter and radiation dominate during the `matter' and `radiation' eras respectively i.e. \Omega_m =1 when w_eff=0 and \Omega_r=1 when w_eff=1/3. We have found that the cosmological dynamical system constrained to obey the LCDM cosmic history has four critical points in each era which correspondingly lead to four forms of f(R). One of them is the usual general relativistic form f(R)=R-2\Lambda. The other three forms in each era, reproduce the LCDM cosmic history but they do not satisfy requirement b. stated above. Only one of these forms (different from general relativity) is found to be an attractor of the dynamical cosmological evolution. It has (\Omega_DE=1, \Omega_r=0, \Omega_m=0) throughout the evolution. Its phase space trajectory is numerically obtained.


I. INTRODUCTION
There is accumulating observational evidence based mainly on Type Ia supernovae standard candles [1] and also on standard rulers [2,3] that the universe has entered a phase of accelerating expansion at a recent cosmological timescale. This expansion implies the existence of a repulsive factor on cosmological scales which counterbalances the attractive gravitational properties of matter on these scales. There have been several theoretical approaches [4,5] towards the understanding of the origin of this factor. The simplest such approach assumes the existence of a positive cosmological constant which is small enough to have started dominating the universe at recent times. The predicted cosmic expansion history in this case (assuming flatness) is H(z) 2 = ȧ a 2 = H 2 0 Ω 0m (1 + z) 3 + Ω 0r (1 + z) 4 + Ω Λ (1.1) where Ω 0r = ρ rad ρcrit ≃ 10 −4 is the present energy density of radiation normalized over the critical density for flatness ρ crit . Also Ω 0m = ρm ρcrit ≃ 0.3 is the normalized present matter density and Ω Λ = 1 − Ω 0m − Ω 0r is the normalized energy density due to the cosmological constant. This model provides an excellent fit to the cosmological observational data [2] and has the additional bonus of simplicity and a single free parameter. Despite its simplicity and good fit to the data this model fails to explain why the cosmological constant is so unnaturally small as to come to dominate the universe at recent cosmological times. This fine tuning problem is known as the coincidence problem.
In an effort to address this problem two classes of models have been proposed: The first class assumes that general relativity is a valid theory on cosmological scales and attributes the accelerating expansion to a dark energy component which has repulsive gravitational properties due to its negative pressure. The role of dark en-ergy is usually played by a minimally coupled to gravity scalar field called quintessence [6]. Alternatively, the role of dark energy can be played by various perfect fluids (eg Chaplygin gas [7]), topological defects [8], holographic dark energy [9] etc. The second class of models attributes the accelerating expansion to a modification of general relativity on cosmological scales which converts gravity to a repulsive interaction at late times and on cosmological scales. Examples of this class of models include scalar-tensor theories [10,11], f (R) modified gravity theories [12], braneworld models [13] etc. An advantage of models in this class is that they naturally allow [10,14] for a superaccelerating expansion of the universe where the effective dark energy equation of state w = p ρ crosses the phantom divide line w = −1. Such a crossing is consistent with current cosmological data [15].
Most of the models in both classes require the existence of arbitrary new degrees of freedom whose role is usually played by effective scalar fields. This is not a welcome feature because the degrees of freedom are to some extend arbitrary with respect to either their origin and/or their dynamical properties. Their predictive power is therefore usually dramatically diminished.
A partial exception to this rule is provided by modified f (R) theories of gravity. In these theories the Ricci scalar R in the general relativistic Lagrangian is replaced by an arbitrary function f (R) leading to an action of the form where L m and L rad are the Lagrangian densities of matter and radiation and we have set 8πG = 1. These f (R) theories arise in a wide range of different frameworks: In quantum field theories in curved spacetime [16], in the low energy limit of the D = 10 superstring theory [17], in the vacuum action for the Grand Unified Theories (GUTs) etc.
It has been demonstrated [18] that for appropriate forms of f (R) the action (1.2) can naturally produce accelerating expansion at late times in accordance with SnIa data [19]. The advantage of these theories is that no extra arbitrary degree of freedom is introduced and the accelerating expansion is produced by the Ricci scalar (dark gravity) whose physical origin is well understood. On the other hand, the main disadvantage of these theories is that (like most modified gravity theories) they are seriously constrained by local gravity experiments [20,21,22]. In fact it can be shown [20] that f (R) models are equivalent to scalar-tensor theories with vanishing Brans Dicke parameter (ω = 0) and a special type of potential. Since solar system tests of general relativity imply ω > 4 × 10 4 [23], these theories can only be consistent with observations if they are associated with a large (infinite) effective mass of the scalar R. It has been shown [24] that specific forms of the function f (R) can provide an infinite effective mass needed to satisfy solar system constraints and can also produce late time accelerating expansion.
The reduction of f (R) theories to a special class of scalar-tensor theories implies that in principle the reconstruction of f (R) from a particular cosmic history H(z) can be performed in a similar way as in the case of scalartensor theories [10,14]. However, the non-existence of a Brans Dicke parameter requires some modifications of the reconstruction methods especially when the reconstruction extends through the whole cosmic history through the radiation and matter eras. The dynamical systems approach followed in the present study illustrates these modifications The construction of cosmologically viable models incorporating late accelerating expansion based on f (R) theories has been an issue of interesting debate during the past year. This debate originated from Ref. [25] which demonstrated that f (R) theories that behave as a power of R at large or small R are not cosmologically viable because they have the wrong expansion rate during the matter dominated era (a ∼ t 1/2 instead of a ∼ t 2/3 ). This conclusion was challenged in Ref. [26] which claimed that wide classes of f (R) gravity models including matter and acceleration phases can be phenomenologically reconstructed by means of observational data. The debate continued with the recent Ref. [27] where a detailed and general dynamical analysis of the cosmological evolution of f (R) theories was performed. It was shown that even though most functional forms of f (R) are not cosmologically viable due to the absence of the conventional matter era required by data, there are special forms of f (R) that can be viable (consistent with data) for appropriate initial conditions.
In the present study we perform a generic model independent analysis of f (R) theories. Instead of specifying various forms of f (R) and finding the corresponding cosmological dynamics, we specify the cosmological dynamics to that of the ΛCDM cosmology and search for a possible corresponding form of f (R). We thus attempt to reconstruct f (R) from the background cosmological dynamics. In particular we consider the general autonomous system for cosmological dynamics of f (R) theories and study the dynamics of f (R) using as input a ΛCDM cosmic expansion history. Our study is performed both analytically (using the critical points and their stability) and numerically by explicitly solving the dynamical system. The results of the two approaches are in good agreement since the numerical evolution of f (R) follows the evolution of the 'attractor' (stable critical point) of the system for most initial conditions. As we point out in the next section however the physical significance of this 'attractor' should be interpreted with care since it is an artifact of the allowed perturbations in the form of the physical law f (R).
The structure of this paper is the following: In the next section we derive the autonomous system for the cosmological dynamics of f (R) theories. Using as input a particular cosmic history H(z) (eg ΛCDM) we show how can this system be transformed so that its solution provides the dynamics and functional form of f (R). We also study the dynamics of this transformed system analytically by deriving its critical points and their stability during the three eras of the cosmic background history (radiation, matter and deSitter). We find that there are 'attractor' critical points for each era which allow an analytical prediction of the dynamics of the system. We also confirm this analytical prediction by a numerical solution of the system demonstrating that the evolution of the system is independent of the initial conditions. In section III we use the solution of the above system to reconstruct the cosmological evolution and functional form of the function f (R). We also demonstrate the agreement between the analytical and numerical reconstruction of f (R). Finally in section IV we conclude, summarize and refer to future prospects of this work.

II. DYNAMICS OF f (R) COSMOLOGIES
We consider the action (1.2) describing the dynamics of f (R) theories in the Jordan frame [12]. In the context of flat Friedman-Robertson-Walker (FRW) universes the metric is homogeneous and isotropic ie and variation of the action (1.2) leads to the following dynamical equations which are the generalized Friedman equations where F ≡ df dR and ρ m , ρ rad represent the matter and radiation energy densities which are conserved according toρ In order to study the cosmological dynamics implied by equations (2.2), (2.3) we express them as an autonomous system [27] of first order differential equations. To achieve this, we first write (2.2) in dimensionless form as We now define the dimensionless variables x 1 , ..., x 4 as 10) where in (2.10) we have used the fact that and we can associate x 4 with Ω r and x 1 + x 2 + x 3 ≡ Ω DE with curvature dark energy (dark gravity). Defining also Ω m ≡ ρm 3F H 2 we can write equation (2.6) as We may now use (2.7) to express (2.3) as Also, differentiating x 4 of (2.11) with respect to N we have where we have made use of (2.5). Similarly, differentiating (2.9) with respect to N we find 19) and , R implies derivative with respect to R. Finally differentiating (2.10) with respect to N we find The autonomous dynamical system (2.15), (2.18), (2.20), (2.17) is the general dynamical system that describes the cosmological dynamics of f (R) theories. It has been extensively studied in Ref. [27] for various cases of f (R) (or equivalently various forms of m) and was found to lead to a dynamical evolution that in most cases is incompatible with observations since it involves no proper matter era. Some forms of f (R) however were found to lead to a cosmological evolution that is potentially consistent with observations. In order to investigate such cases in more detail we follow a different approach. Instead of investigating the above autonomous system for various different behaviors of m(f (R)) we eliminate m from the system by assuming a particular form for H(N ) (ie x 3 (N ) (see (2.10))) consistent with cosmological observations. Once x 3 (N ) is known we can solve (2.20) for x1x3 m and substituting in (2.18) we find which along with (2.15) and (2.17) consist a new dynamical system which is independent of m. The study of this system will be our focus in what follows.
The results of our analysis do not rely on the use of any particular form of x 3 (N ) (ie H(z)). They only require that the universe goes through the radiation era (high redshifts), matter era (intermediate redshifts) and acceleration era (low redshifts). The corresponding total effective equation of state is For the sake of definiteness however, we will assume a specific form for H(z) corresponding to a ΛCDM cosmology (1.1) which in terms of N takes the form We can use (2.10) and (2.24) to find x 3 (N ) as 25) The crucial generic properties of x 3 (N ) are its values at the radiation, matter and deSitter eras: (2.28) where N rm ≃ −ln Ω0m Ω0r and N mΛ ≃ − 1 3 ln ΩΛ Ω0m are the N values for the radiation-matter and matter-deSitter transitions. For Ω 0m = 0.3, Ω 0r = 10 −4 we have N rm ≃ −8, N mΛ ≃ −0.3. The transition between these eras is model dependent but rapid and it will not play an important role in our analysis.
It is straightforward to study the dynamics of the system (2.15), (2.21), (2.17) by setting x ′ i = 0 to find the critical points and their stability in each one of the three eras corresponding to (2.26)-(2.28). Notice that even though this dynamical system is not autonomous at all times it can be approximated as such during the radiation, matter and deSitter eras when x 3 is approximately constant. The critical points and their stability are shown in Table I.
The stability analysis of Table I assumes that x 3 = const and therefore it is not identical to the full stability analysis where x 3 would be allowed to vary. The usual stability analysis of cosmological dynamical systems assumes a particular cosmological model (eg a form of f (R) or m) and in the context of this 'physical law', the stability of cosmic histories H(N ) is investigated. In this context clearly a stable cosmic history is the one preferred by the model.
In the reconstruction approach however the stability analysis has a very different meaning. Here we do not fix the model f (R) ('physical law'). Here we fix the cosmic history and allow the physical law f (R) to vary in order to predict the required cosmic history. Thus our stability analysis concerns the 'physical law' f (R) and not the particular cosmic history. Since the physical law is usually fixed by Nature the instabilities we find are not physically relevant but they are only useful to understand analytically the phase space trajectories we obtain numerically. The physically interesting quantities are the values of the critical points we find in each era in the context of the ΛCDM cosmic history. These tell us the possible physical laws f (R) that can reproduce a ΛCDM cosmic history. As shown in Table I, one of these laws is clearly the general relativistic f (R) = R − 2Λ.
The important point to observe in Table I is that in each era there are four critical points one of which correspond to the general relativistic f (R) = R − 2Λ. Three of the four critical points in each era are not stable. This however does not imply that these points are not cosmologically relevant. These instabilities are not instabilities of the trajectory H(N ) (which we keep fixed) but of the forms of f (R) which is allowed to vary. Thus they are not so relevant physically since in a physical context f (R) is assumed to be fixed a priori. The 'attractor' critical points of Table I are relevant only for technical reasons since they allow a comparison between a numerical evolution and an analytical prediction of the evolution of the system. In a more realistic situation where the perturbations of the 'physical law' f (R) would be turned off, all critical points would correspond to valid reconstructions of a ΛCDM cosmic history.
If we allow for f (R) perturbations (but not of H(N ) perturbations), the evolution of the system is determined by just following the evolution of the 'attractors' of Table  I through the three eras. This evolution is presented in Table II showing the 'attractors' in each era. We stress however that this is not necessarily a preferred cosmological trajectory for the reasons described above. The 'standard' critical points are shown in Table III and they are the only critical points that have in addition to the correct expansion rate properties, the required values of (Ω r , Ω m , Ω DE ) in each era. As shown in the next section these saddle critical points reconstruct the general relativistic f (R) ie f (R) = R − 2Λ. It is therefore clear that nonlinear f (R) theories can produce an observationally acceptable cosmic history but not with the required values of (Ω r , Ω m , Ω DE ) in each era. We should stress that our analysis has not excluded the possibility of physical values of (Ω 0r , Ω 0m , Ω DE ) in the case of cosmic histories oscillating around the anticipated w ef f in each era or a w ef f that is continuously evolving. These special cases however maybe severely constrained observationally.
To confirm the dynamical evolution implied by the 'attractors' of Table I, we have performed a numerical analysis of the dynamical system (2.15), (2.21), (2.17) using the ansatz (2.25) for x 3 with Ω 0m = 0.3 and Ω 0r = 10 −4 . This ansatz for x 3 (N ) leads to the w ef f (N ) shown in Fig.  1. We have set up the system initially, close to the 'standard' radiation era saddle point (0, 0, 0, 1) and allowed it to evolve. As seen in Fig. 2 and Fig. 3 much before the onset of the matter era (N ≡ N rr ≃ −25 < −8 ≃ N rm ) the slow (but non-zero) evolution of x 3 (N ) forces the phase space trajectory to depart from the saddle point (0, 0, 0, 1) and head towards the radiation era stable 'attractor' (−4, 5, 0, 0) where it stays throughout the rest of the radiation era (w ef f ≃ 1 3 ).         Table III is also shown (continuous green line). The points A1, A2, A3 correspond to the 'attractors' of each era (radiation, matter and deSitter respectively) while the points B1, B2, B3 correspond to the 'standard' critical points of each era (see Tables II and III). Notice that on the projection of Fig. 3b the 'attractor' points A2, A3 coincide with the 'standard' critical points B2, B3.
Subsequently, when x 3 (N ) enters the matter era (w ef f = 0) at N rm ≃ −8, the trajectory follows the evolution of the 'attractor' fixed point and heads towards the matter era 'attractor' (−3.386, 3.886, 0.5, 0) ignoring the saddle point (0, −1/2, 1/2, 0) of the 'standard' matter era. Finally when the matter era is over, the trajectory heads towards the deSitter 'attractor' (−1, 0, 2, 0) which is also distinct from the 'standard' deSitter saddle point (0, −1, 2, 0). Notice that the deSitter 'attractor' is inconsistent with observations due to the implied large variation of the effective Newton's constant G(N ) = 1 F (N ) even though this inconsistency could be ameliorated by 'chameleon' type mechanisms [29]. The evolution of (Ω r , Ω m , Ω DE ) corresponding to the phase space trajectories of Figs. 2 and 3 is shown in Fig. 4a. Notice that Ω m = 0 throughout the 'attractor' evolution of the system and the w ef f = 0 of the matter era is induced by curvature dark gravity excitations.
We have also tested initial conditions in the matter era starting the evolution on the saddle point (0, − 1 2 , 1 2 , 0) corresponding to the 'standard' matter era. In this case we also ignore radiation setting Ω 0r = 0. We get an evolution of the system (see Figs 2,3,4b) which stays on the 'standard' matter era Ω m = 1 for about 3 expansion times but before the onset of the acceleration era it gets absorbed by the 'attractor' towards the nonstandard deSitter critical point (−1, 0, 2, 0).
The above evolution along the 'attractor' critical points is a result of the 'physical law' f (R) perturbations. We can also reproduce trajectories that go through critical points that are not stable by turning off these perturbations. For example we can recover the saddle critical point sequence (0, 0, 0, 1) → (0, by fixing x 1 = 0 in the system (2.15), (2.21), (2.17) and reducing it to the system which can be easily solved using the ansatz (2.25) (see Fig. 2 and Fig. 3). The alternative approach of solving the decoupled pair (2.18), (2.17) does not lead to the correct result because the f (R) perturbations are not turned off and the constraint is not respected in this case.
It is straightforward to reconstruct the functions f (R) that correspond to the saddle general relativistic trajectory (2.29) and to the 'attractor' sequence of Table II. The functional forms of f (R) may also be reconstructed on any one of the critical points of Table I. These tasks are undertaken in the next section.

III. RECONSTRUCTION OF f (R)
We now reconstruct the form of the function f (R) that corresponds to each one of the critical points of the system shown in Table I. This reconstruction is effectively an approximation of f (R) in the neighborhood of each critical point. It is particularly useful because most of the dynamical evolution takes place close to the fixed points. Consider a critical point of the form (x 1 ,x 2 ,x 3 ,x 4 ). Using (2.8) we find where F 0 is a constant. We may eliminate N in favor of R using the input form of H(N ) (equation (2.24) in (2.12) to obtain (setting H 2 0 = 1) which leads to and by integration we get where f 0 is an integration constant. Expressing (3.4) in terms of N using (3.2) we obtain It is now straightforward to use the expressions for f (N ) R(N ) and H(N ) to find x 2 (N ) (equation (2.9)), x 3 (N ) (equation (2.25)) and x 4 (N ) (equation (2.11)). We thus find 3.6) and while x 3 (N ) is given by (2.25). Using equations (3.6), (2.25) and (3.7) we may verify thex 2 ,x 3 ,x 4 values of each critical point by considering the appropriate range of N in each era and the corresponding value ofx 1 . By demanding consistency with the values of Table I we may obtain the values of the constants f 0 and F 0 .
As an example let's consider the sequence (2.29) corresponding to the 'standard' cosmological eras (Table III). It is easy to see, usingx 1 = 0 and the appropriate range of N in (3.6), (2.25) and (3.7) that we obtain the correct values forx 2 ,x 3 ,x 4 in the radiation and matter eras for any value of F 0 , f 0 . In the deSitter era (N >> 1) , Ω rad (dotted) and ΩDE (short-dashed) parameters for 'standard' radiation era initial conditions (Fig4a) and 'standard' matter era initial conditions (Fig4b). In both cases the 'attractor' critical points of Table II rapidly take over and drag the trajectories.
the value of f 0 needs to be fixed to get agreement with x 2 = −1 of Table I. In particular from (3.6) we find Using now (3.9) and settingx 1 = 0 in (3.4) we reconstruct the expected result which is valid for all three eras since in this sequence the value ofx 1 remains constant. In a similar way we may reconstruct f (R) for any critical point in one of the three eras. We have therefore extended previous studies showing that f (R) theories can only be viable in very restricted cases by showing that even these restricted cases can not reproduce a viable ΛCDM cosmology where w ef f is constant during the matter and radiation eras and Ω 0r , Ω m take their cosmologically anticipated values. It therefore becomes clear that if the accelerating expansion of the universe is due to physics in the gravitational sector it may probably have to be a more general theory than f (R) modified gravity. Such a theory could very well be scalartensor gravity (or equivalently coupled dark energy [30]) whose cosmological dynamical properties and constraints need to investigated in detail.
In the case of sequence of transitions among critical points which involve different values ofx 1 the reconstruction can be done by either numerical determination of x 1 (N ) or by approximating it as a sequence of step functions. For example, the steps involved in the reconstruction of the 'attractor' trajectory shown in Fig. 2 The numerical solution x 1 (N ) of Fig. 2a can be approximated as a piecewise constant function with values determined by the corresponding 'attractors' of each cosmological era and by the initial conditions ie 3.15) (where N rr ≃ −25) thus leading to an analytical approximation for F (N ). The resulting form of ln(F (N )) in both the numerical reconstruction and its analytical approximation is shown in Fig. 5 (5a: 'standard' radiation era initial condition, 5b: 'standard' matter era initial condition).
2. Use equation (2.9) to find f (N ) from F (N ) ie where H(N ) is given by (2.24), x 2 (N ) is numerically obtained and shown in Fig. 2b and F (N ) is obtained in the previous step. As in the case of x 1 (N ), x 2 (N ) can be analytically approximated as 3.20) using the corresponding 'attractors' to obtain an analytical expression for f (N ).  N )) in the numerical reconstruction (dashed lines) and its analytical approximation using both the interpolating steps 1-3 (dotted lines) and the application of the analytical expression (3.5) valid in each era (thick green lines). The agreement between the three approaches is very good. 5a: 'standard' radiation era initial condition, 5b: 'standard' matter era initial condition. The form of f (R) in the numerical reconstruction (dashed lines) and its analytical approximation using both the interpolating steps 1-3 (dotted lines) and the application of the analytical expression (3.4) valid in each era (continuous green lines). 6a: 'standard' radiation era initial condition with the continuous green lines corresponding to radiation era (larger R) and matter era (smaller R). The deSitter era is not shown since it corresponds to a single point (constant R). 6b: 'standard' matter era initial condition with the continuous green line corresponding to matter era. The deSitter era is not shown since it corresponds to a point (constant R).
3. The resulting form of f (N ) can then be combined with equation (3.2) for R(N ) to reconstruct the function f (R). The resulting form of f (R) is shown in Fig. 6 for both the numerical reconstruction and its analytical approximation (6a: radiation era initial conditions, 6b: matter era initial conditions).
We can fit the reconstructed f (R) of Figs. 6a and 6b to the analytic form of equation (3.4) for each era respectively so as to find the parameters F 0 ,x 1 and f 0 . The results are shown in Table IV. Notice that the best fit values ofx 1 coincide with the corresponding 'attractor' critical points of Table II as expected. This verifies the validity of the reconstructed f (R) expression from (3.4). A similar reconstruction analysis can be performed for any other sequence of critical points. As discussed in section II any such sequence is equally interesting cosmologically since the existence of the 'attractor' is an artifact of the f (R) perturbations.

IV. CONCLUSION-OUTLOOK
We have shown analytically and numerically that nonlinear f (R) gravity theories can reproduce the background expansion history H(z) indicated by observations even when f (R) does not reduce to general relativity at early times. In that case the universe gets dominated by dark gravity during its evolution as opposed to radiation or matter. This result relies on the values of all the critical points we found assuming only that the radiation era corresponds to a constant effective equation of state parameter w ef f = 1 3 while for the matter era we have w ef f = 0.
Our analysis indicates f (R) models can be viable if f (R) deviates from general relativity at early times. Thus a viable f (R) theory must satisfy one of the fol-lowing: • Either f (R) reduces to general relativity at early times, but departs from general relativity at late times (a well known case [12]).
• Or dark gravity in the forms derived in our paper mimics radiation or matter at both the background level and the perturbative level. The later would clearly require a separate analysis of perturbations of the model.
Numerical Analysis: The mathematica files with the numerical analysis of this study may be found at http://leandros.physics.uoi.gr/frlcdm/frlcdm.htm or may be sent by e-mail upon request.